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Abstract 


The  objective  of  this  workshop  was  to  cross-fertilize  new  developments 
in  acoustics  and  signal  processing  in  the  field  of  ocean  acoustics  and 
tomography,  laboratory  ultrasonics,  seismology  and  medical  ultrasonics. 
The  workshop  was  held  in  Cargese,  Corsica  (France).  Groups  of  15 
leading  scientific  researcheers  presented  lectures  in  each  of  their 
respective  fields. 


Research  Summary 


The  workshop  combined  a  tutorial  approach  in  providing  a  common 
interdisciplinary  basis  with  presentations  and  discussions  of  the  latest 
research  results.  Contributions  (13)  of  each  of  the  discussions  were 
compiled  in  a  workshop  proceedings  entitled,  “Imaging  of  Complex 
Media  with  Acoustic  and  Seismic  Waves”  and  submitted  to  Springer- 
Verlag  GmBH  &  Co.,  Germany  for  publishing. 

A  chapter  entitled,  “Ocean  Acoustics,  Matched  Field  Processing  and 
Phase  Conjugation”  by  W.  A.  Kuperman  and  D.  R.  Jackson  is  included 
from  the  book  as  submittal  of  the  reporting  requirements  under  this  grant. 
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INTRODUCTION 

This  chapter  treats  ocean  acoustics  and  various  applications  of  signal  processing,  phase  conjugation  and 
tomography  to  ocean  acoustics.  Phase  conjugation  is  related  to  the  standard  principles  of  passive  signal 
processing  and  to  the  more  recent  demonstrations  of  time  reversal  mirrors  in  the  ocean.  Here  we  will 
overview  relevant  physics  of  ocean  acoustics,  present  a  passive  signal  processing  method  called  matched 
field  processing  and  then  discuss  the  related  concept  and  implementation  of  time  reversal  methods  in 
ocean  acoustics. 

1  Review  of  Ocean  Acoustics 

The  ocean  is  a  waveguide,  bounded  above  by  a  pressure  release  surface  and  below  by  a  visco-elastic 
medium.  The  physical  oceanographic  parameters,  as  ultimately  represented  by  the  ocean  sound  speed 
structure,  make  up  the  index  of  refraction  of  the  water  column  waveguide.  The  combination  of  water 
column  and  bottom  properties  leads  to  a  set  of  generic  sound  propagation  paths  descriptive  of  most 
propagation  phenomena  in  the  ocean. 
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Figure  1 :  Generic  sound  speed  profiles 


1.1  QUALITATIVE  DESCRIPTION  OF  OCEAN  SOUND  PROPAGATION 

First  we  summarize  aspects  of  oceanography  that  impact  propagation  paths  before  we  go  into  more 
detail  on  the  propagation  itself.  Figure  1  illustrates  a  typical  set  of  sound  speed  profiles  indicating 
greatest  variability  near  the  surface  as  a  function  of  season  and  time  of  day.  In  a  warmer  season  (or 
warmer  part  of  the  day),  the  temperature  increases  near  the  surface  and  hence  the  sound  speed  decreases 
with  depth.  In  nonpolar  regions,  the  oceanographic  properties  of  the  water  near  the  surface  result  from 
mixing  activity  originating  from  the  air-sea  interface.  This  near  surface  mixed  layer  has  a  constant 
temperature  (except  in  calm,  warm  surface  conditions  as  described  above).  In  this  isothermal  mixed 
layer,  the  sound  speed  profile  can  increase  with  depth  due  to  the  pressure  gradient  effect.  This  is  the 
“surface  duct”  region. 

1.1.1  Sound  Propagation  Paths  in  the  Ocean 

Below  the  mixed  layer  is  the  thermocline  in  which  the  temperature  decreases  with  depth  and  the  sound 
speed  decreases  with  depth.  Below  the  thermocline,  the  temperature  is  constant  and  the  sound  speed 
increases  because  of  increasing  pressure.  Therefore,  there  exists  a  depth  between  the  deep  isothermal 
region  and  the  mixed  layer  with  a  minimum  in  sound  speed;  this  depth  is  often  referred  to  as  the  axis 
of  the  deep  sound  channel. 

However,  in  polar  regions,  the  water  is  coldest  near  the  surface  and  hence  the  minimum  sound  speed 
is  at  the  ocean/air(or  ice)  interface  as  indicated  in  Fig.  1  In  continental  shelf  regions  (shallow  water) 
with  water  depths  on  the  order  of  a  few  hundred  meters,  only  the  upper  region  of  the  sound  speed 
profile  in  Fig.  1,  which  is  dependent  on  season  and  time  of  day,  affects  sound  propagation  in  the  water 
column. 

Figure  2  is  a  contour  display  of  the  sound  speed  structure  of  the  North  and  South  Atlantic  [1] 
with  the  axis  of  the  deep  sound  channel  indicated  by  the  heavy  dashed  line.  Note  the  geographic(and 
climatic)  variability  of  the  upper  ocean  sound  speed  structure  and  the  stability  of  this  structure  in  the 
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Figure  2:  Sound  speed  contours  at  5  m/s  intervals  taken  from  the  North  and  South  Atlantic  along  30.50 
degrees  West.  Dashed  line  indicates  axis  of  deep  sound  channel  (from  Northrup  and  Colburn,  1974). 


deep  isothermal  layer.  For  example,  as  explained  above,  the  axis  of  the  deep  sound  channel  becomes 
shallower  toward  both  poles,  eventually  going  to  the  surface.  Figure  3  is  a  schematic  of  the  basic  types 
of  propagation  in  the  ocean  resulting  from  the  sound  speed  profiles  (indicated  by  the  dashed  lines) 
discussed  in  the  last  section.  These  sound  paths  can  be  understood  from  a  simplified  statement  of 
Snell’s  Law:  sound  bends  locally  toward  regions  of  low  sound  speed  (or  sound  is  “trapped”  in  regions 
of  low  sound  speed).  Paths  A  and  B  correspond  to  surface  duct  propagation  where  the  minimum  sound 
speed  is  at  the  ocean  surface  (or  at  the  bottom  of  the  ice  cover  for  the  Arctic  case).  Path  C,  depicted  by 
a  ray  leaving  a  deeper  source  at  a  shallow  horizontal  angle  propagates  in  the  deep  sound  channel  whose 
axis  is  at  the  shown  sound  speed  minimum.  This  local  minimum  tends  to  become  more  shallow  towards 
polar  latitudes  converging  to  the  Arctic  surface  minimum.  Hence  for  mid  latitudes,  sound  in  the  deep 
channel  can  propagate  long  distances  without  interacting  with  lossy  boundaries;  propagation  via  this 
path  has  been  observed  over  distances  of  thousands  of  kilometers.  Also,  from  the  above  description 
of  the  geographical  variation  of  the  acoustic  environment  combined  with  Snell’s  Law,  we  can  expect 
that  shallow  sources  coupling  into  the  water  column  at  polar  latitudes  will  tend  to  propagate  more 
horizontally  around  an  axis  which  becomes  deeper  toward  the  mid  latitudes.  Path  D,  which  is  at 
slightly  steeper  angles  than  those  associated  with  path  C,  is  convergence  zone  propagation,  a  spatially 
periodic  (~  35  -  65  km)  refocusing  phenomenon  producing  zones  of  high  intensity  near  the  surface  due 
to  the  upward  refracting  nature  of  the  deep  sound  speed  profile.  Referring  back  to  Fig.  1,  there  may 
be  a  depth  in  the  deep  isothermal  layer  at  which  the  sound  speed  is  the  same  as  it  is  at  the  surface. 
This  depth  is  called  the  “critical  depth”  and,  in  effect,  is  the  lower  limit  of  the  deep  sound  channel. 
A  “positive  critical  depth”  specifies  that  the  environment  supports  long  distance  propagation  without 
bottom  interaction,  whereas  ‘hiegative”  implies  that  the  bottom  ocean  boundary  is  the  lower  boundary 
of  the  deep  sound  channel.  The  bottom  bounce  path,  E,  which  interacts  with  the  ocean  bottom  is 
also  a  periodic  phenomenon  but  with  a  shorter  cycle  distance  and  a  shorter  total  propagation  distance 
because  of  losses  when  sound  is  reflected  from  the  ocean  bottom.  Finally,  the  right  hand  side  of  Fig.  3 
depicts  propagation  in  a  shallow  water  region  such  as  a  continental  shelf  Here  sound  is  channeled  in  a 
waveguide  bounded  above  by  the  ocean  surface  and  below  by  the  ocean  bottom. 

The  modeling  of  sound  propagation  in  the  ocean  is  further  complicated  because  the  environment 
varies  laterally  (range  dependence),  and  all  environmental  effects  on  sound  propagation  are  dependent 
on  acoustic  frequency  in  a  rather  complicated  way  which  often  makes  the  ray  type  schematic  of  Fig. 
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Figure  3:  Schematic  representation  of  various  types  of  sound  propagation  in  the  ocean. 

3  misleading,  particularly  at  low  frequencies.  Finally,  a  quantitative  understanding  of  acoustic  loss 
mechanisms  in  the  ocean  is  required  for  modeling  sound  propagation.  These  losses  are,  aside  from 
geometric  spreading:  volume  attenuation,  bottom  loss  (i.e.,  a  smooth  water/bottom  interface  is  not  a 
perfect  reflector),  and  surface,  volume  (including  fish)  and  bottom  scattering  loss.  Here  we  will  only 
review  those  aspects  of  bottom  loss  which  impact  propagation  structure. 

1.1.2  Bottom  Loss 

Ocean  bottom  sediments  are  often  modeled  as  fluids  since  the  rigidity  (and  hence  the  shear  speed)  of 
the  sediment  is  usually  considerably  less  than  that  of  a  solid,  such  as  rock.  In  the  latter  case,  which 
applies  to  the  “ocean  basement”  or  the  case  where  there  is  no  sediment  overlying  the  basement,  the 
medium  must  be  modeled  as  an  elastic  solid  which  means  it  supports  both  compressional  and  shear 
waves. 

Reflectivity,  the  amplitude  ratio  of  reflected  and  incident  plane  waves  at  an  interface  separating  t\vo 
media,  is  an  important  measure  of  the  effect  of  the  bottom  on  sound  propagation.  For  an  interface 
between  two  fluid  semi-infinite  halfspaces  with  density  pi  and  sound  speed  c,-,  z  =  1,2,  as  shown  in 
Fig.  4(a)  (assuming  a  harmonic  time  dependence  of  exp(— iwi)),  the  reflectivity  is  given  by 


p2ku^plk2z 

(1) 

with 

hiz  =  (w/ci)  sin  Qi  =  ki  sin  ;  i  =  1 , 2. 

(2) 

The  incident  and  transmitted  grazing  angles  are  related  by  Snell’s  Law, 

kjL  =  cc>s  6i  =  k2  cos  ^2, 

(3) 

where  the  incident  grazing  angle  Bi  is  also  equal  to  the  angle,  of  the  reflected  plane  wave.  Jl{0)  is 
also  referred  to  as  the  Rayleigh  Reflection  Coefficient  and  has  unit  magnitude  (total  internal  reflection) 
when  the  numerator  and  denominator  of  Eq.  1  are  complex  conjugates.  This  occurs  when  k2z  is  purely 


Kuperman  and  Jackson:  Ocean  Acoustics 


5 


a)  -V 

exp  [I  (kj,  •  r  +  kizZ)]  R  (0)exp  (i  (k^  •  r  -  k^jZ)) 


Figure  4:  The  reflection  and  transmission  process.  Grazing  angles  are  defined  relative  to  the  horizontal. 

imaginary  and,  using  Snell’s  Law  for  determining  62  in  terms  of  the  incident  grazing  angle,  we  obtain 
the  critical  grazing  angle  below  which  there  is  perfect  reflection, 

cos^c=ci/c2,  (4) 

so  that  a  critical  angle  can  exist  only  when  the  speed  in  the  second  medium  is  higher  than  that  of  the 
first. 

Using  Eq.  2,  Eq.  1  can  be  rewritten  as 

_  P2C2/sin^2  -Pici/sin^i  ^  Z2  -  Zi 
P2C2/ sm02  d- piCi/ sinOi  ^2  +  ^1’ 

where  the  right  hand  side  is  in  the  form  of  impedances^  Zi{6i)  =  piCi/  sin  ,  which  are  the  ratios  of  the 
pressure  to  the  vertical  particle  velocity  at  the  interface  in  the  medium.  Written  in  this  form,  the 
reflectivity  for  more  complicated  boundaries  follows  in  a  straightforward  manner  [2]. 

In  lossy  media,  attenuation  can  be  included  in  the  reflectivity  formula  by  taking  the  sound  speed  as 
complex  so  that  the  wavenumbers  are  subsequently  also  complex,  ki  ^  ki  +  a,-. 

Figure  4(b)  depicts  a  simple  bottom  loss  curve  derived  from  the  Rayleigh  reflection  coefficient 
formula  where  both  the  densities  and  sound  speed  of  the  second  medium  are  larger  than  those  in  the 
first  medium  with  unit  reflectivity  indicating  perfect  reflection.  For  loss  in  dB,  0  dB  is  perfect  reflecting, 
6  dB  loss  is  an  amplitude  factor  of  one-half,  12  dB  loss  is  one-fourth,  etc.  For  a  lossless  bottom,  severe 
loss  occurs  above  the  critical  angle  in  the  water  column  due  to  transmission  into  the  bottom.  For  the 
lossy  (more  realistic)  bottoms,  only  partial  reflection  occurs  at  all  angles.  With  paths  involving  many 
bottom  bounces  (shallow  water  propagation),  bottom  losses  as  small  as  a  few-tenths  of  a  dB  per  bounce 
accumulate  and  become  significant  because  the  propagation  path  may  involve  many  tens  of  bounces. 
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Figure  5:  Ocean  waveguide  propagation,  (a)  Long  distance  propagation  occurs  within  a  cone  of  2^c- 
(b)  Same  as  (a)  but  with  nonisovelocity  water  column  showing  refraction. 

Path  E  in  Fig.  3,  the  bottom  bounce  path,  often  involves  paths  which  correspond  to  angles  near  or 
above  the  critical  angle;  therefore,  after  a  few  bounces,  the  sound  level  is  highly  attenuated.  On  the 
other  hand,  for  shallow  angles,  many  bounces  are  possible.  Hence,  in  shallow  water,  path  F,  most  of 
the  energy  that  propagates  is  close  to  the  horizontal  and  this  type  of  propagation  is  most  analogous 
to  waveguide  propagation.  In  fact,  as  shown  in  Fig.  5,  there  exists  a  small  cone  from  which  energy 
propagates  long  distances  (Oc  is  typically  10°  ~  20°).  Energy  outside  the  cone  is  referred  to  as  the  near 
field  (or  continuous  spectrum),  which  eventually  escapes  the  waveguide.  The  trapped  field  originating 
from  within  the  cone  is  referred  to  as  the  normal  mode  field  (or  discrete  spectrum)  because  there  are 
a  set  of  angles  corresponding  to  discrete  paths  which  constructively  interfere  and  make  up  the  normal 
(natural)  modes  of  the  shallow  water  environment. 

1.2  SOUND  PROPAGATION  MODELS 

Sound  propagation  in  the  ocean  is  mathematically  described  by  the  wave  equation,  whose  parameters 
and  boundary  conditions  are  descriptive  of  the  ocean  environment.  There  are  essentially  four  types  of 
models  (computer  solutions  to  the  wave  equation)  to  describe  sound  propagation  in  the  sea;  ray  theory, 
the  spectral  method  or  fast  field  program  (FFP),  normal  mode  (NM)  and  parabolic  equation  (PE)  [2]. 
All  of  these  models  allow  for  the  fact  that  the  ocean  environment  varies  with  depth.  A  model  that  also 
takes  into  account  horizontal  variations  in  the  environment  (i.e.,  sloping  bottom  or  spatially  variable 
oceanography)  is  termed  range  dependent.  For  high  frequencies  (a  few  kilohertz  or  above),  ray  theory 
is  the  most  practical.  The  other  three  model  types  are  more  applicable  and  usable  at  lower  frequencies 
(below  a  kilohertz).  Here,  we  will  confine  the  discussion  to  the  normal  mode  model  including  its  mild 
adiabatic  range-dependent  extension  and  the  parabolic  equation  (see  Appendix  A). 
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1*2.1  The  Wave  Equation  and  Boundary  Conditions 

The  wave  equation  for  pressure  in  cylindrical  coordinates  with  the  range  coordinates  denoted  by  r  = 
(x,y)  and  the  depth  coordinate  denoted  by  z  (taken  positive  downward)  for  a  source  free  region  of 
constant  density  p  is 

V 


where  c  is  the  sound  speed  in  the  wave  propagating  medium.  The  wave  equation  is  most  often  solved  in 
the  frequency  domain;  that  is,  a  frequency  dependence  of  exp{—iu3i)  is  assumed  to  obtain  the  Helmholtz 
equation  {K  =  w/c), 

V^p{v,z)’]- K^p{r,z)  =  0.  (7) 

The  most  common  plane  interface  boundary  conditions  encountered  in  underwater  acoustics  are 
described  below:  For  the  ocean  surface  there  is  the  pressure  release  condition  where  the  pressure 
vanishes, 

p=0.  (8) 

The  interface  between  the  water  column  (layer  1)  and  an  ocean  bottom  sediment  (layer  2)  is  often 
characterized  as  a  fluid-fluid  interface.  The  continuity  of  pressure  and  vertical  particle  velocity  at  the 
interface  yields  the  following  boundary  conditions  in  terms  of  pressure: 

1  dpi  _  1  dpi 

T:~d7-7,~dl- 

These  boundary  conditions  applied  to  the  plane  wave  fields  in  Fig.  4(a)  yield  the  Rayleigh  reflection 
coefficient  given  by  Eq.  1. 

The  Helmholtz  equation  for  an  acoustic  field  from  a  point  source  with  angular  frequency  u  is 

V^G{i,z)  +  K\T,z)G(v,z)  =  -6^iT-T,)8{z-z,y,  K^{t,z)= 

where  the  subscript  “s”  denotes  the  source  coordinates.  The  range  dependent  environment  manifests 
itself  as  the  coefficient  z)  of  the  partial  differential  equation  for  the  appropriate  sound  speed  profile. 

The  range  dependent  bottom  type  and  topography  appear  as  boundary  conditions.  The  acoustic  field 
from  a  point  source,  G(r),  is  either  obtained  by  solving  the  boundary  value  problem  of  Eq.  10  (spectral 
method  or  normal  modes)  or  by  approximating  Eq.  10  by  an  initial  value  problem  (ray  theory,  parabolic 
equation) . 

1.2.2  Normal  Mode  Model  (NM) 

For  a  range-independent  environment,  we  assume  a  pressure  field  solution  of  Eq.  10  of  the  form 


1 

G(r,  z)  =  —  /  <Pk  g(k,  z,  z,)  exp[ik-(r  -  r,)], 
J -oo 


which  then  leads  to  the  equation  for  the  depth  dependent  Green’s  function,  g(k,  z,z,), 


(11) 
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We  can  solve  Eq.  12  for  g  using  a  normal  mode  expansion  of  the  form 

(13) 

where  the  quantities  Un  are  eigenfunctions  of  the  following  eigenvalue  problem: 

* 

-j^+[K'‘{z)-kl]un{z)  =  0.  (14) 

The  eigenfunctions,  Un,  are  zero  at  =  0,  satisfy  the  local  boundary  conditions  descriptive  of  the  ocean 
bottom  properties  and  satisfy  a  radiation  condition  for  ^  oo.  They  form  an  orthonormal  set  in  a 
Hilbert  space  with  weighting  function  p{z\  the  local  density.  The  range  of  discrete  eigenvalues  of  Eq.  14 
is  given  by  the  condition 

min[/L(2:)]  <kn<  max[/f(2;)].  (15) 

These  discrete  eigenvalues  correspond  to  discrete  angles  within  the  critical  angle  cone  in  Fig.  5  such 
that  specific  waves  constructively  interfere.  The  mode  functions  form  a  complete  set  (for  simplicity  we 
omit  discussion  of  the  continuous  spectrum  though  a  good  approximation  is  to  use  a  set  of  discrete 
mode  functions  obtained  from  a  waveguide  extended  in  depth  and  terminated  by  a  pressure  release  or 
rigid  boundary) 

W„(i)K„(z,)  ,,  , 

2^  — ^771 —  =  0{Z-Zs),  (16) 

all  modes  rV 

and  satisfy  the  orthonormality  condition 


Um(z)Un(z) 


dz  —  ^nm ) 


where  6nm  is  the  Kronecker  delta  symbol.  The  eigenvalues  kn  typically  also  have  a  small  imaginary 
part  On ,  which  serves  as  the  modal  attenuation  representative  of  all  the  losses  in  the  ocean  environment 
(see  [2]  for  the  formulation  of  normal  mode  attenuation  coefficients). 

Solving  Eq.  10  using  the  normal  mode  expansion  given  by  Eq.  13  yields  (for  the  source  at  the  origin). 


G(r,  z)—  ^  (18) 

The  asymptotic  form  of  the  Hankel  function  can  be  used  in  the  above  equation  to  obtain  the  well  known 
normal  mode  representation  of  a  cylindrical  (axis  is  depth)  waveguide: 


G(r,  z)  = 


p{zs){STTry/^ 


exp(-i7r/4)  ^ 


■Un{Zs)Un{z) 


exp{iknr). 


Equation  19  is  a  far  held  solution  of  the  wave  equation  and  neglects  the  continuous  spectrum  (k^  < 
min[/^(2:)]  of  Ineq.  15)  of  modes.  For  purposes  of  illustrating  the  various  portions  of  the  acoustic  field, 
we  note  that  kn  is  a  horizontal  wave  number  so  that  a  ^^ray  angle”  associated  with  a  mode  with  respect 
to  the  horizontal  can  be  taken  to  be  ^  =  cos^^[kn/ K{z)].  For  a  simple  waveguide  the  maximum  sound 
speed  is  the  bottom  sound  speed  corresponding  to  min[/C(2:)].  At  this  value  of  K{z)y  we  have  from 
Snell’s  Law  0  =  0^^  the  bottom  critical  angle.  In  effect,  if  we  look  at  a  ray  picture  of  the  modes,  the 
continuous  portion  of  the  mode  spectrum  corresponds  to  rays  with  grazing  angles  greater  than  the 
bottom  critical  angle  of  Fig.  4(b)  and  therefore  outside  the  cone  of  Fig.  5.  This  portion  undergoes 
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Figure  6:  Ray  and  normal  mode  theory,  (a)  Sound  speed  profile  (b)  Ray  trace  (c)  Normal  modes  (d) 
Propagation  calculations. 


severe  loss.  Hence,  we  note  that  the  continuous  spectrum  is  the  near  (vertical)  field  and  the  discrete 
spectrum  is  the  (more  horizontal,  profile  dependent)  far  field  (falling  within  the  cone  in  Fig.  5). 

The  advantages  of  the  NM  procedure  are  that:  the  solution  is  available  for  all  source  and  receiver 
configurations  once  the  eigenvalue  problem  is  solved;  it  is  easily  extended  to  moderately  range-dependent 
conditions  using  the  adiabatic  approximation;  it  can  be  applied  (with  more  effort)  to  extremely  range- 
dependent  environments  using  coupled  mode  theory.  However,  it  does  not  include  a  full  representation 
of  the  near  field. 


1.2.3  Adiabatic  Mode  Theory 


All  of  the  range  independent  normal  mode  “machinery”  developed  for  environmental  ocean  acoustic 
modeling  applications  can  be  adapted  to  mildly  range  dependent  conditions  using  adiabatic  mode  theory. 
The  underlying  assumption  is  that  individual  propagating  normal  modes  adapt  (but  do  not  scatter  or 
“couple”  into  each  other)  to  the  local  environment.  The  coefficients  of  the  mode  expansion,  a„  in  Eq.  13, 
now  become  mild  functions  of  range,  i.e.,  a„(k)  — +  a„(k,r).  This  modifies  Eq.  19  as  follows: 


G(r,z)  = 


p{zs)(SiTry/^ 


exp(-ix/4)  ^  exp(i*„r). 


where  the  range-averaged  wave  number  (eigenvalue)  is 


and  the  kn(r^)  are  obtained  at  each  range  segment  from  the  eigenvalue  problem  Eq.  14  evaluated  at  the 
environment  for  that  particular  range  along  the  path.  The  quantities  and  t;„  are  the  sets  of  modes 
at  the  source  and  the  field  positions,  respectively. 
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Figure  7:  Relationship  between  FFP,  NM  and  PE  computations,  (a)  FFP  Green’s  function  from  Eq.  12. 
(b)  NM,  FFP  and  PE  propagation  results  showing  some  agreement  in  near  field  and  complete  agreement 
in  far  field. 

Simply  stated,  the  adiabatic  mode  theory  leads  to  a  description  of  sound  propagation  such  that 
the  acoustic  field  is  a  function  of  the  modal  structure  at  both  the  source  and  the  receiver  and  some 
average  propagation  conditions  between  the  two.  Thus,  for  example,  when  sound  emanates  from  a 
shallow  region  where  only  two  discrete  modes  exist  and  propagates  into  a  deeper  region  with  the  same 
bottom  (same  critical  angle),  the  two  modes  from  the  shallow  region  adapt  to  the  form  of  the  first  two 
modes  in  the  deep  region.  However,  the  deep  region  can  support  many  more  modes;  intuitively,  we 
therefore  expect  the  resulting  two  modes  in  the  deep  region  will  take  up  a  smaller  more  horizontal  part 
of  the  cone  of  Fig.  5  than  they  take  up  in  the  shallow  region.  This  means  that  sound  rays  going  from 
shallow  to  deep  tend  to  become  more  horizontal  which  is  consistent  with  a  ray  picture  of  downslope 
propagation.  Finally,  fully  coupled  mode  theory  for  range  dependent  environments  has  been  developed 
[3]  but  requires  extremely  intensive  computation. 

1.3  QUANTITATIVE  DESCRIPTION  OF  PROPAGATION 

All  of  the  models  described  above  attempt  to  describe  reality  and  to  solve  in  one  way  or  another  the 
Helmholtz  equation.  They  therefore  should  be  consistent  and  there  is  much  insight  to  be  gained  from 
understanding  this  consistency.  The  models  ultimately  compute  propagation  loss  which  is  taken  as  the 
decibel  ratio  (see  Appendix  B)  of  the  pressure  at  the  field  point  to  a  reference  pressure,  typically  one 
meter  from  the  source. 

Figure  6  shows  convergence  zone  type  propagation  for  a  simplified  profile.  The  ray  trace  in  Fig.  6(b) 
shows  the  cyclic  focusing  discussed  in  Section  1.1.  The  same  profile  is  used  to  calculate  the  normal 
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Figure  8:  Model  and  data  comparison  for  a  range  dependent  case,  (a)  Profiles  and  ray  trace  for  a  case 
of  a  surface  duct  disappearing,  (b)  250  Hz  PE  and  2  kHz  Ray  trace  comparisons  with  data. 

modes  shown  in  Fig.  6(c)  which  when  summed  according  to  Eq.  19  yield  the  same  cyclic  pattern  as  the 
ray  picture.  Figure  6(d)  shows  both  the  normal  mode  (wave  theory)  and  ray  theory  result.  Ray  theory 
exhibits  sharply  bounded  shadow  regions  as  expected  whereas  the  normal  mode  theory,  which  includes 
diffraction,  shows  that  the  acoustic  field  does  exist  in  the  shadow  regions  and  the  convergence  zones 
have  structure. 

Normal  mode  models  sum  the  discrete  modes  which  roughly  correspond  to  angles  of  propagation 
within  the  cone  of  Fig.  5.  The  spectral  (FFP)  method  can  include  the  full  field,  discrete  plus  continuous, 
the  latter  corresponding  to  larger  angles.  The  spectral  model  involves  solving  Eq.  12  numerically  for  a 
continuum  of  wavenumbers,  k,  without  the  normal  mode  expression  of  Eq.  13  (see  Fig.  7(a)).  Therefore, 
the  consistency  we  expect  between  the  normal  mode  and  the  spectral  method  and  the  physics  of  Fig.  5 
is  that  the  continuous  portion  of  the  spectral  solution  decays  rapidly  with  range  so  that  there  should 
be  complete  agreement  at  long  ranges  between  normal  mode  and  spectral  solutions.  The  Lloyd’s  mirror 
effect,  a  near  field  effect,  should  also  be  exhibited  in  the  spectral  solution  but  not  the  normal  mode 
solution.  These  aspects  are  apparent  in  Fig.  7(b).  The  PE  solution  appears  in  Fig.  7(b)  and  is  in  good 
agreement  with  the  other  solutions  but  with  some  phase  error  associated  with  the  average  wavenumber 
that  must  be  chosen  in  the  split  step  method.  The  PE  solution,  which  contains  part  of  the  continuous 
spectrum  including  the  Lloyd  mirror  beams,  is  more  accurate  than  the  normal  mode  solution  at  short 
range;  more  recent  PE  results  [27]  can  be  made  arbitrarily  accurate  in  the  forward  direction. 

Range  dependent  results  [2]  are  shown  in  Fig.  8.  A  ray  trace,  a  ray  trace  field  result,  a  PE  result  and 
data  are  plotted  together  for  a  range  dependent  sound  speed  profile  environment.  The  models  agree 
with  the  data  in  general  with  the  exception  that  the  ray  results  predict  too  sharp  a  leading  edge  of  the 
convergence  zone. 

Upslope  propagation  is  modeled  with  the  PE  in  Fig.  9.  As  the  field  propagates  upslope,  sound  is 
dumped  into  the  bottom  in  what  appear  to  be  discrete  beams  [4].  The  fiat  region  has  three  modes 


Kuperman  and  Jackson:  Ocean  Acoustics  12 


Figure  9:  Relation  between  upslope  propagation  (from  PE  calculation)  showing  individual  mode  cutoff 
and  energy  dumping  in  the  bottom,  and  a  corresponding  ray  schematic. 

and  each  is  cut-off  successively  as  sound  propagates  into  shallower  water.  The  ray  picture  also  has  a 
consistent  explanation  of  this  phenomenon.  The  rays  for  each  mode  become  steeper  as  they  propagate 
upslope.  When  the  ray  angle  exceeds  the  critical  angle  the  sound  is  significantly  transmitted  into  the 
bottom.  The  locations  where  this  takes  place  for  each  of  the  modes  is  identified  by  the  three  arrows. 

2  Matched  Field  Processing 

Spatial  sampling  of  a  sound  field  is  usually  done  by  an  array  of  transducers,  although  the  synthetic 
aperture  array,  in  which  a  sensor  is  moved  through  space  to  obtain  measurements  in  both  the  time  and 
space  domains,  is  an  important  exception.  Spatial  sampling  is  analogous  to  temporal  sampling  with 
the  sampling  interval  replaced  by  the  sensor  spacing  vectors.  The  Nyquist  criterion  requires  that  the 
sensor  spacing  be  at  least  twice  the  spatial  wavenumber  of  the  measured  sound  field.  Aliasing  and  signal 
reconstruction  follow  the  same  principles  as  for  temporal  data  except  that  ambiguity  can  arise  in  the 
reconstructed  wavenumber  components  due  to  symmetries  in  the  array  geometry.  The  simplest  example 
of  array  processing  is  phase  shading  in  the  frequency  domain  (or  time  delay  in  the  time  domain)  to 
search  in  bearing  for  a  plane  wave  signal  and  is  referred  to  as  plane  wave  beamforming,  or  delay  and 
sum  beamforming  in  the  time  domain.  This  form  of  beamforming  is  discussed  elsewhere  and  we  will 
move  on  to  a  process  referred  to  as  matched  field  processing. 

Matched  field  processing  (MFP)  [5]  is  the  three-dimensional  generalization  of  the  conventional  lower- 
dimensional  plane  wave  beamformer.  The  one  or  two  dimensions  in  the  latter  case  are  bearing  and/or 
elevation  and  the  matching  is  done  to  plane  waves.  Localization  in  this  case  refers  to  determining  a 
direction.  The  generalized  matched  field  beamformer  matches  the  measured  field  at  the  array  with 
replicas  of  the  expected  field  for  all  source  locations.  These  replicas  are  derived  from  propagation 
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Figure  10:  Matched  field  processing,  (a)  schematic  of  the  procedure,  (b)  block  diagram  of  the  processor. 

models  as  discussed  above  or  in  more  detail  in  Ref.  [2].  The  unique  spatial  structure  of  the  field  permits 
localization  in  range,  depth  and  azimuth  depending  on  the  array  geometry  and  complexity  of  the  ocean 
environment.  The  complex  interference  pattern  of  the  acoustic  field  is  a  function  of  the  source  location 
and  this  pattern  can  be  “matched,”  In  terms  of  rays  we  can  say  that  the  refractive  properties  of 
the  waveguide  generate  a  pattern  of  arrival  angles  that  can  also  be  matched.  The  process  is  shown 
schematically  in  Fig.  10.  The  process  consists  of  systematically  placing  a  test  point  source  at  each 
point  of  a  search  grid,  computing  the  acoustic  field  (replicas)  at  all  the  elements  of  the  array  and  then 
correlating  this  modeled  field  with  the  data  from  the  real  point  source  whose  location  is  unknown.  When 
the  test  point  source  is  collocated  with  the  true  point  source,  the  correlation  will  be  a  maximum.  The 
two  main  factors  which  limit  performance  of  MFP  are  noise,  natural  and  man  made,  and  the  ability  to 
accurately  model  the  ocean  acoustic  environment. 

Now  construct  the  replica  field,  w(a)  on  the  array  for  each  candidate  position,  a  from  one  of  the 
numerical  solutions  to  the  acoustic  wave  equation,  Eq,  10.  The  replica  vector  is  normalized,  i.e.,  it  has 
unit  length,  so  that  the  output  of  the  processor  is  the  source  strength, 

”<•>=1511'  <“> 

where  G(a)  is  the  vector  whose  zth  element  is  the  solution  to  the  waveguide  equation,  e.g.,  Eq.  19,  at 
the  field  (hydrophone)  position  (r,j?)  =  (ra^JSj)  for  a  source  at  a. 

Rather  than  find  just  the  direction  of  the  source  from  the  array,  we  can  search  for  its  actual  location 
by  matching  the  received  data,  d(atrue)  from  the  true  source  location  against  solutions  of  the  wave 
equation.  The  matching  process  is  then  similar  to  plane  wave  beamforming  with  the  plane  waves  being 
replaced  by  solutions  of  the  wave  equation.  The  output  of  this  matched  field  process,  denoted  5(a),  at 
each  point  in  space  a  is  given  by 

S(a)  =  |w^(a)d(a(r«e)p.  (23) 

This  equation  can  be  rewritten  in  terms  of  the  cross-spectral  density  matrix  of  the  array  data,  K(airue)  = 

d(ajrue)d^  (^true)i 

S(a)  =  w^(a)K(a,rue)w(a).  (24) 
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Figure  11:  Matched  field  processing,  (a)  Environment  and  geometry  for  shallow  water  example,  (b) 
Matched-field  ambiguity  surfaces  for  a  130  dB  source  (r5=10  km,  2^=50  m).  Bartlett  (top)  and  MVDP 
(bottom). 


There  are  an  assortment  of  interesting  issues  concerning  the  actual  measurement  of  cross-spectral  density 
matrix,  particularly  when  adaptive  processing  is  utilized.  These  issues  are  addressed  in  the  chapter  by 
B  agger  oer. 

The  peak  of  the  output  of  the  beamformer,  S(a),  is  at  atme^  S(a)  is  also  referred  to  as  the  ambiguity 
function  (or  surface)  of  the  matched  held  processor  because  it  also  contains  ambiguous  peaks  which  are 
analogous  to  the  sidelobes  of  a  conventional  plane  wave  beamformer.  Sidelobe  suppression  can  often  be 
accomplished  by  adaptive  processing  which  is  also  discussed  in  the  Baggeroer  chapter.  Here  we  write 
down  the  form  of  one  of  the  adaptive  processors:  the  minimum  variance  distortionless  processor 

SMVDp(a)  =  [w^(a)K“^(ajr„e)w(a)]“\  (25) 

where  again  we  note  from  above  that  there  are  measurement  issues  associated  with  estimating  an 
invertible  cross  spectral  density  matrix. 

A  vertical  receive  array  simulation  example  of  the  Bartlett  and  MVDP  MFP  processors  for  an  ocean 
acoustic  waveguide  with  a  high  signal  to  noise  ratio  is  shown  in  Fig.  11.  In  MFP,  mismatch  between 
replicas,  that  are  derived  from  a  usually  incomplete  knowledge  of  the  environment,  and  observed  data  is 
a  significant  problem  that  has  received  much  attention.  Matched  field  tomography  (MFT)searches  for 
the  environmental  parameters  controlling  the  propagation  (for  example  the  index  of  refraction  which 
may  be  a  spatially  dependent  coefficient  of  the  wave  equation)  rather  than  source  location.  As  with  the 
case  of  plane  wave  beamforming,  there  are  an  assortment  of  processors  that  can  be  used  for  MFP/MFT 

[5] 

There  have  been  two  approaches  to  deal  with  environmental  mismatch  in  MFP.  The  first  has  been 
to  constructing  processors  which  are  tolerant  to  environmental  uncertainty  [6,  7,  8].  These  processors 
trade  resolution  for  lack  of  precise  knowledge  of  the  environment.  The  second  approach  is  focalization 
[9]  which  uses  nonlinear  optimization  methods  to  simultaneously  search  for  the  environment  and  the 
source  position.  All  of  these  methods  have  limitations,  but  under  various  conditions  they  have  been 
successfully  applied  to  experimental  data. 
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3  Phase  Conjugation  in  the  Ocean 

Phase  conjugation  (PC)  [10]  and  its  time  domain  equivalent  referred  to  as  a  time  reversal  mirror  (TRM) 
is  a  process  that  has  been  first  demonstrated  in  nonlinear  optics.  The  time  reversal  mirror  in  ultrasonics 
[11,  12]  is  the  subject  of  Ch.  ???.  Aspects  of  phase  conjugation  as  applied  to  underwater  acoustics  also 
have  been  explored  recently  [13,  14,  15,  16].  The  Fourier  conjugate  of  phase  conjugation  is‘time  reversal; 
implementation  of  such  a  process  over  a  finite  spatial  aperture  results  in  a  “time  reversal  mirror.”  This 
section  will  begin  with  a  description  of  the  basic  theory  and  implementation  in  the  ocean  based  on  the 
ocean  acoustics  discussion  above.  This  will  be  followed  by  a  discussion  of  some  of  the  more  complex 
issues  associated  with  the  real  ocean  such  as  fluctuations,  etc.  Some  of  the  material  given  below  will 
overlap  and  supplement  Ch.  ???. 

3.1  Basic  Properties  of  Phase  Conjugation 

Phase  conjugation  takes  advantage  of  reciprocity  which  is  a  property  of  wave  propagation  in  a  static 
medium  and  is  a  consequence  of  the  invariance  of  the  linear  lossless  wave  equation  to  time  reversal.  In 
the  frequency  domain,  time  reversal  corresponds  to  conjugation  invariance  of  the  Helmholtz  equation. 
The  property  of  reciprocity  allows  one  to  retransmit  a  time  reversed  version  of  a  multipath  dispersed 
probe  pulse  back  to  its  origin,  arriving  there  time  reversed,  with  the  multipath  structure  having  been 
undone  [17,  18].  This  process  is  equivalent  to  using  the  ocean  as  a  matched  filter  since  the  probe  pulse 
arrival  has  embedded  in  it  the  transfer  function  of  the  medium.  This  process  can  be  extended  further 
by  receiving  and  retransmitting  the  probe  signal  with  a  source-receive  array.  Depending  on  the  spatial 
extent  of  the  array,  the  above  process  results  in  some  degree  of  spatial  focusing  of  the  signal  at  the 
origin  of  the  probe  signal. 

A  time  reversal  mirror  (TRM)  can  therefore  be  realized  with  a  source-receiver  array.  The  incident  sig¬ 
nal  is  received,  time  reversed  and  transmitted  from  sources  contiguous  with  the  receiving  hydrophones. 
The  time  reversal  can  be  accomplished  in  a  straightforward  way,  for  example,  by  using  the  rewind  out¬ 
put  of  an  analog  tape  recorder  or  by  a  simple  program  that  reverses  a  digitized  segment  of  a  received 
signal. 

Phase  conjugation  or  the  implementation  of  a  TRM  in  the  ocean  will  be  shown  to  be  related  to 
matched  field  processing  (MFP)  which  requires  detailed  knowledge  of  the  environment.  Phase  con¬ 
jugation  is  an  environmentally  self-adaptive  process  which  may  therefore  have  significant  applications 
to  localization  and  communications  in  complicated  ocean  environments.  Though  the  “effective”  ocean 
environment  must  remain  static  over  the  turn  around  time  of  the  PC  process,  ocean  variability  on  time 
scales  shorter  than  the  turn  around  time  might  be  compensated  for  with  feedback  algorithms.  However, 
an  understanding  of  relevant  ocean  time  scales  vis  a  vis  the  stability  of  the  PC  process  will  be  required. 

As  noted  earlier,  phase  conjugation  in  the  ocean  shares  many  features  in  common  with  ultrasonic 
time  reversal.  These  general  aspects  will  be  discussed  first  and  then  those  aspects  that  demand  treat¬ 
ment  of  the  ocean  waveguide  will  be  covered.  The  generic  phase  conjugation  focusing  experiment  is 
depicted  in  Fig.  12.  A  point  acoustic  source,  denoted  the  “probe  source”  transmits  a  signal  which  is 
received  by  the  elements  of  an  array  capable  of  both  reception  and  transmission  (the  “source-receiver 
array”,  or  “SRA”).  The  signal  received  by  each  element  is  time-reversed  and  retransmitted  on  the 
same  element  (or  perhaps  on  a  nearby  source  element  if  the  receiving  elements  are  unsuited  as  acoustic 
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Figure  12:  Schematic  of  the  generic  phase  conjugation  focusing  experiment. 

sources).  The  held  produced  by  the  array  approximates  a  time-reversed  version  of  the  probe  source 
field,  consequently  it  focuses  at  the  location  of  the  probe  source. 

3.1.1  Idealized  Arrays 

In  order  to  assess  the  degradation  in  focusing  due  to  environmental  effects,  it  is  helpful  to  have  a  reference 
point  for  comparison.  Thus,  we  will  first  consider  ideal  arrays  capable  of  providing  the  sharpest  possible 
focus.  While  one  might  think  of  the  ideal  case  as  producing  a  simple  time-reversed  version  of  the  field 
due  to  the  probe  source,  this  ideal  is  unreachable,  even  in  principle.  To  see  this,  consider  a  point 
probe  source,  which  produces  a  field  that  approaches  infinity  as  one  approaches  the  source  location. 
This  infinite  field  does  not  satisfy  the  source-free  wave  equation  at  the  source  location,  therefore  it  is 
impossible  to  produce  such  a  field  using  an  array  whose  elements  are  at  a  distance  from  the  source 
location. 

These  considerations  lead  to  the  question,  how  sharp  can  the  focus  of  the  time-reversed  field  be?  It 
is  useful  to  first  consider  the  ideal  case,  in  which  a  continuous  SRA  surrounds  a  probe  source,  as  shown 
in  Fig.  13. 

The  probe  source,  situated  at  the  location  Xps  emits  a  time-harmonic  field  of  angular  frequency  w. 
We  will  take  the  source  to  have  unit  strength,  obeying  the  inhomogeneous  Helmholtz  equation 

V^Ga,(r,  rp,)  +  A:(r)^Ga,(r ,  Tp,)  =  -(5(r  -  rp,)  .  (26) 

Note  that  this  Green’s  function  incorporates  all  the  complexity  of  the  problem  at  hand.  That  is,  it 
satisfies  appropriate  boundary  conditions  such  as  those  at  the  sea  surface  and  bottom  and  includes  all 
the  effects  of  inhomogeneity  in  the  sound  speed  in  the  position-dependent  wavenumber  A:(r)  =  w/c(r). 
At  a  general  point  r'  on  its  surface,  the  ideal  SRA  receives  the  field  Gu;(r',rp,)  of  the  probe  source, 
conjugates  it,  and  retransmits  to  a  point  r  in  the  interior  according  to  Huygen’s  principle  as  embodied 
in  the  Helmholtz-Kirchhoff  integral  [30].  The  resulting  field  at  r  is 
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Figure  13:  Schematic  of  ideal,  closed  source-receiver  array  (SRA). 


Pp,(r;a;)  =  ^  [G„(r,  rOG:(r',rp,)  +  VG„(r,r')  •  VG:(r^rp,)]  dS'  .  (27) 

As  Porter  [30]  has  shown,  in  the  lossless  case  (that  is,  when  k(r)  is  real),  the  field  produced  by  the 
ideal  array  can  be  expressed  exactly  as 


Ppc(r;u;)  =  22^m{Ga,(r,rp,)}  .  (28) 

By  taking  the  imaginary  part  of  Eq.  26  it  can  be  seen  that  Ppc(r)  rpai  w)  is  a  solution  of  the  sourceless 
(homogeneous)  Helmholtz  equation.  Thus,  it  is  seen  that  the  ideal  SRA  produces  afield  which  is  similar 
to  the  original  field,  but  is  source-free  in  the  focus  region,  as  it  must  be.  The  most  remarkable  aspect  of 
this  ideal  case  is  the  focus  field  is  independent  of  the  shape  and  location  of  the  closed  array,  as  long  as 
it  surrounds  the  original  source.  We  will  see  that  this  ideal  behavior  has  implications  in  more  practical 
situations. 

How  good  is  the  ideal  focus?  In  the  case  where  the  medium  is  homogeneous  in  the  vicinity  of  the 
original  source  and  in  which  backscattering  due  to  boundaries  and  inhomogeneities  produces  a  negligible 
field  near  the  original  source,  the  original  field  can  be  approximated  by  the  free-space  Green’s  function, 


^  e®p(j^|r-rp,|) 

G.(r,rp,)«  •  (29) 

While  this  field  is  singular  at  r  =  Vpg ,  its  imaginary  part  is  not,  as  shown  in  Fig  14.  It  can  be  seen  from 
the  figure  that  the  ideal  focus  has  a  width  comparable  to  the  wavelength.  This  ideal  is  not  generally 
realized  in  more  practical  situations.  Note  that  the  focus  field  is  oscillatory.  It  is,  in  fact,  a  standing 
wave  created  by  the  interference  of  waves  launched  from  all  directions,  since  the  ideal  SRA  surrounds 
the  focus. 


3.1.2  The  open  array 
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Figure  14:  Comparison  of  Source  and  Ideal  Focus  Fields 

In  more  realistic  cases,  the  array  will  not  be  closed  and  will  not  surround  the  original  source.  The 
ideal  array  provides  insight  into  such  cases,  as  one  can,  as  a  thought  experiment,  divide  the  closed 
SRA  into  two  parts,  as  in  Fig.  15,  the  smaller  of  these  parts  representing  an  actual  two-dimensional 
SEA.  Consider  deforming  this  array,  while  leaving  the  other  portion  undeformed.  We  know  from  the 
previous  result  that  the  focus  field  will  not  change  as  a  result  of  this  deformation.  The  focus  field  in  this 
case  is  the  sum  of  the  fields  from  the  two  parts,  and,  as  one  part  of  the  array  is  fixed,  its  contribution 
to  the  field  must  be  fixed  as  well.  It  follows  that  the  field  of  the  smaller  part  must  be  unchanging, 
too.  Thus,  an  ideal  open  two-dimensional  SRA  should  display  an  interesting  invariance  property:  Its 
focus  field  will  be  independent  of  deformations  of  the  array  provided  the  outer  edge  of  the  array  is  held 
fixed.  Of  course,  this  array  is  rather  special;  it  is  continuous  and  employs  both  the  incident  field  and  its 
normal  derivative.  Nevertheless,  we  shall  see  that  this  invariance  survives  at  least  as  an  approximation 
in  realistic  situations. 

It  is  possible  to  extract  an  additional  useful  piece  of  information  from  this  example.  Consider  an 
infinite,  planar  SRA  in  a  homogeneous  acoustic  medium.  It  can  be  considered  as  one-half  of  a  closed 
array  formed  from  two  planes  placed  symmetrically  with  respect  to  the  probe  source.  In  the  plane 
parallel  to  the  arrays  and  passing  through  the  source,  symmetry  dictates  that  each  array  produces  the 
same  field,  thus  the  infinite  planar  SRA  must  produce  a  sin(kR)/kR  field  in  the  plane  through  the 
source,  where  R  is  the  distance  from  the  source. 

3.1.3  The  discrete  array 

In  practice,  creating  an  array  to  approximate  even  an  open  version  of  the  array  represented  by  Eq.  13  is 
very  difficult,  as  this  array  must  be  capable  of  sensing  both  the  incident  field  and  its  normal  derivative. 
Furthermore,  it  must  include  both  monopole  sources  (to  produce  G(r,r',a;)),  and  dipole  sources  (to 
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Figure  15:  Demonstration  of  shape  invariance  of  ideal  open  arrays. 

produce  VG(r,  r',a;)).  The  general  problem  of  using  arrays  with  various  combinations  of  monopole  and 
dipole  sources  is  discussed  by  Cassereau  and  Fink  [31].  In  low-frequency  underwater  acoustics,  arrays 
are  usually  composed  of  discrete  elements  that  are  omnidirectional  or  nearly  so.  Such  arrays  are  more 
accurately  represented  by  a  set  of  discrete  monopole  receivers/sources  as  illustrated  in  Fig.  12. 

In  the  harmonic  case,  such  an  array  produces  the  following  field  in  response  to  a  point  probe  source 
of  unit  strength. 

J 

Fpc(r;u>)  =  ^G2,u,(r,rj)Gl^(rj-,rp,)  .  (30) 

j=i 

The  sum  is  over  the  J  elements  of  the  SRA  whose  position  vectors  are  denoted  ry .  Propagation  from  the 
probe  source  to  the  SRA  is  described  by  the  Green’s  function  Gi^w(ry,rp5),  while  propagation  from  the 
SRA  to  the  field  point  is  described  by  G2,a;(r,  r„).  The  subscripts  1  and  2  allow  for  the  possibility  that 
time  variation  of  the  ocean  might  cause  changes  in  the  Green’s  function  between  the  probe  and  phase- 
conjugate  transmission  cycles.  During  either  propagation  cycle,  the  ocean  is  assumed  to  be  ‘frozen”  in 
the  sense  that  it  behaves  as  a  time-invariant  linear  system.  In  this  view,  the  Green’s  function  is  the 
frequency-dependent  system  transfer  function  for  acoustic  propagation  between  any  two  points  in  the 
ocean.  It  instructive  to  consider  the  focus  field  at  the  location  of  the  probe  source.  Suppressing  the 
subscripts  1  and  2  in  Eq.  30,  the  field  at  r  =  Xpg  is 

J 

P pcij^ps  \  |Gaj(rpg ,  ry)|^  .  (31) 

i=i 

In  Eq.  31,  reciprocity  has  been  used  in  the  form  Ga;(rp5,ry)  =  Ga,(ry,rp,),  which  neglects  gradients  in 
density.  The  field  at  the  focus  is  a  sum  of  positive  terms,  that  is,  the  transmissions  from  each  element 
of  the  SRA  arrive  exactly  in  phase  at  the  probe  source  location.  Even  though  each  term  in  this  sum 
is  likely  to  have  complicated  dependence  on  frequency,  u;,  owing  to  modal  interference  (equivalently, 
owing  to  multipath  propagation),  the  sum  will  tend  to  have  less  severe  frequency  dependence  than  the 
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Figure  16:  Experimental  setup  of  phase  conjugation  experiment. 

individual  terms.  This  is  because  each  element  of  the  SRA  is  subject  to  different  modal  interference, 
with  destructive  interference  tending  to  occur  at  different  frequencies  for  each  element.  This  lack  of 
frequency  dependence  at  the  source  location  implies  that  a  time-reversed  pulse  will  suffer  little  multipath 
distortion  at  the  focus.  See  Sec.  3.2.2  for  a  discussion  of  this  matter  in  the  time  domain.  Note  that 
the  field  at  the  focus  is  proportional  to  the  sum  over  the  SRA  of  the  squared  probe  source  pressure  at 
each  element.  If  the  SRA  samples  the  nearly  the  entire  water  depth,  this  sum  will  tend  to  be  stable 
and  independent  of  details  of  propagation  in  the  intervening  medium,  as  the  sum  is  a  measure  of  the 
total  power  propagating  through  the  waveguide.  Thus,  if  one  imagines  perturbing  the  environment 
by  altering  the  seafloor  roughness  or  by  introducing  slight  random  changes  in  the  sound  speed  of  the 
medium,  one  would  not  expect  the  field  value  at  the  focus  to  change  appreciably.  This  is  another 
indication  of  the  robust  nature  of  phase-conjugate  focusing. 

3.2  Background  Theory  and  Simulation  for  Phase  Conjugation/TRM  in  the 
Ocean 

The  theory  of  phase  conjugation  vis  a  vis  ocean  acoustics  already  has  been  presented  [19,  20].  Here  we 
briefly  review  salient  issues  using  the  basic  geometry  of  an  ocean  TRM  experiment  [19,  21]  as  shown 
schematically  in  Fig.  16.  The  figure  also  indicates  the  types  of  environmental  measurements  that  were 
made.  The  TRM  was  implemented  by  a  77  m  source-receiver  array  (SRA)  in  125  m  deep  Avater  The  SRA 
consisted  of  20  hydrophones  with  20  contiguously  located  sources  with  a  nominal  resonance  frequency 
of  445  Hz.  The  sources  were  operated  at  a  mean  nominal  165  dB  source  level.  The  received  signals 
from  a  probe  source  (PS)  were  digitized,  time  reversed  and  after  being  converted  back  to  analog  form, 
retransmitted.  A  vertical  46  element  receive  array  (VRA)  spanning  90  m  located  initially  6.3  km  from 
the  SRA  at  the  PS  range. 
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Figure  17:  Single  frequency  simulation  of  phase  conjugation  for  the  geometry  of  Fig.  16  for  a  probe 
source  located  at  a  depth  of  40  m  and  a  range  of  6.3  km.  (a)  Sound  speed  profile.  The  density,  p  and 
attenuation,  a  (in  dB/ wavelength)  of  the  bottom  two  layers  are  also  given,  (b)  Simulation  for  a  20 
element  SRA.  Note  the  sharp  focus  in  depth,  (c)  Simulation  for  only  the  bottom  10  elements  of  the 
SRA. 


3.2.1  Harmonic  Point  Source 

For  a  unit  harmonic  source  of  angular  frequency  uj  located  at  (vpsyZps),  the  pressure  field,  Gco{R\  ZjyZps)^ 
at  the  j-th  receiver  element  of  the  SRA  from  the  point  source  PS  in  Fig.  16  is  determined  from  Eq.  19, 


Gu;{r;ZyZps)  = 


p{zp,)(8Trry/^ 


exp(-i7r/4)^ 


^n(gp»)^n(^) 


exp{iknr). 


For  a  vertical  line  of  discrete  sources,  Eq.30  (with  the  subscripts  1  and  2  dropped,  because  time  depen¬ 
dence  of  the  medium  is  neglected),  takes  the  form 

J 

Ppc{r,z-,u)  =  'Y^Gu{r-,z,Zj)Gl,{R]Zj,Zp,),  (33) 

i=i 

where  R  is  the  horizontal  distance  of  the  SRA  from  PS  and  r  is  the  horizontal  distance  from  the  SRA  to 
a  field  point.  We  can  implement  Eq.  33  with  propagation  models  other  than  the  normal  mode  example 
we  are  using.  For  mildly  range-dependent  environments,  we  can  use  the  adiabatic  mode  theory  for 
Guj(r\ZyZj)  described  in  Section  1.2.3.  For  more  highly  range-dependent  environments,  the  parabolic 
equation  model  is  not  only  more  appropriate,  it  is  intuitively  suitable  in  the  sense  that  it  uses  a  range 
marching  algorithm;  hence,  it  is  a  direct  implementation  of  back  propagation.  This  model  is  described 
in  Appendix  A  and  we  note  here  that  it  has  recently  been  applied  to  some  interesting  medical  ultrasonics 
problems  [22]. 

Note  that  the  magnitude  squared  of  the  right  hand  side  (r.h.s.)  of  Eq.  33  is  the  unnormalized 
(see  Eq.  22)  ambiguity  function  of  the  matched  held  processor  of  Eq.  23  where  the  data,  d,  is  given 
by  Gw{R'i  )  ^ps)  and  the  unnormalized  replica  field  by  Gu,{r]ZjZj),  In  effect,  the  process  of  phase 
conjugation  is  an  implementation  of  matched  field  processing  where  the  ocean  itself  is  used  to  construct 
the  replica  field.  Or,  alternatively,  matched  field  processing  simulates  the  experimental  implementation 
of  phase  conjugation  in  which  a  source/receive  array  is  used.  To  demonstrate  that  Ppd'f'y  z)  focuses  at 
the  position  of  the  probe  source,  (i^,  Zp^),  we  simply  substitute  Eq.  32  into  Eq.  33  which  specifies  that 
we  sum  over  all  modes  and  array  sources 


Ppcir,z-,w)FaY,J212 


)‘^n(^j‘  )^n(^p3) 

Pi^j  )/^(^pi)‘N/^m kjiVR 


exp  i{kn 
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For  an  array  which  substantially  spans  the  water  column  and  adequately  samples  most  of  the  modes, 
we  may  approximate  the  sum  of  sources  as  an  integral  and  invoke  orthonormality  as  specified  by  Eq.  17 . 
Then  the  sum  over  j  selects  out  modes  m  —  n  and  Eq.  34  becomes 


Ppc(r,2;w) 

m 


p{^^ps)kYfiy/vR 


expikm{r  —  R)- 


(35) 


The  individual  terms  change  sign  rapidly  with  mode  number.  However,  for  the  field  at  PS,  r  —  R,  the 
closure  relation  of  Eq.  16  can  be  applied  approximately  (we  assume  that  the  A:„’s  are  nearly  constant 
over  the  interval  of  the  contributing  modes)  with  the  result  that  Ppci^y  z)  «  S{z  —  Zps).  Figure  17(b)  is 
a  simulation  of  the  phase  conjugation  process  using  Eq.  33  for  a  probe  source  at  40  m  depth  and  at  a 
range  of  6.3  km  from  a  20  element  SRA  as  specified  in  Fig.  16  verifying  the  above  discussion.  Measured 
range- dependent  bathymetry,  bottom  properties  and  sound  speed  profile  were  used  as  the  input  to  an 
adiabatic  mode  model  (section  1.2.3).  Notice  that  the  focusing  in  the  vertical  is  indicative  of  the  closure 
property  of  the  modes.  As  a  matter  of  fact,  for  an  SRA  with  substantially  fewer  elements,  we  see  that 
the  focusing  still  is  relatively  good.  For  example,  Fig.  17(c)  also  shows  a  result  for  the  bottom  10 
elements  of  the  SRA  which  are  below  the  thermocline. 


3.2.2  Pulse  Excitation 

In  the  actual  experiments,  a  50  ms  pure-tone  pulse  with  center  frequency  445  Hz  was  used  for  the 
probe  transmission.  We  can  Fourier  synthesize  the  above  results  to  examine  phase  conjugation  for 
pulse  excitation.  Here,  in  the  context  of  this  experiment,  we  remind  the  reader  that  phase  conjugation 
in  the  frequency  domain  is  equivalent  to  time  reversal  in  the  time  domain.  The  jth  element  of  the  SRA 
receives  the  following  time-domain  signal,  given  by  Fourier  synthesis  of  the  solution  of  Eq.  28, 

PiR,  zj  ■,t)  =  J  GUR-,  ,  (36) 

where  5(w)  is  the  Fourier  transform  of  the  probe  source  pulse.  This  expression  incorporates  all  waveg¬ 
uide  effects,  including  time  elongation  due  to  multipath  propagation.  For  convenience,  take  the  time 
origin  such  that  P(R,Zj]t)  =  0  outside  the  time  interval  (0,r).  Then  the  time-reversed  signal  that 
will  be  used  to  excite  the  jth  transmitting  element  of  the  SRA  is  P{R^Zj\T  —  t)  such  that  T  >  2r. 
This  condition  is  imposed  by  causality;  the  signal  has  to  be  completely  received  before  it  can  be  time 
reversed.  Then 

P{R,zf,T-t)  =  j G^{R-,Zj,Zp,)S{w)e-'^^'^-^^duj  =  j [Gl{R-, Zj , S' duj  ,  (37) 

where  the  sign  of  the  integration  variable,  w,  has  been  reversed  and  the  conjugate  symmetry  of  the 
frequency-domain  Green’s  function  and  probe  pulse  has  been  used.  The  quantity  in  brackets  in  Eq.  37 
is  the  Fourier  transform  of  the  signal  received  by  the  jth  SRA  receiver  element  after  time  reversal  and 
time  delay.  Hence  there  is  an  equivalence  of  time  reversal  and  phase  conjugation  in  their  respective 
time  and  frequency  domains. 

Noting  that  the  bracketed  quantity  in  Eq.  37  is  the  frequency  domain  representation  of  the  signal 
retransmitted  by  the  jth  element  of  the  SRA,  Fourier  synthesis  can  be  used  to  obtain  the  time-domain 
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Figure  18:  Simulation  of  a  445  Hz,  50  ms  transmitted  pulse  for  the  geometry  in  Fig.  12  for  a  probe 
source  located  at  a  depth  of  40  m.  (a)  Pulse  received  on  the  SR  A  at  range  of  6.3  km  from  PS.  There  is  a 
temporal  dispersion  of  about  75  ms  and  significant  energy  throughout  the  water  column,  (b)  The  focus 
of  the  time  reversed  pulse  at  the  VRA.  There  is  pulse  compression  back  to  the  original  transmitted  50 
ms  duration  as  well  as  spatial  focusing  in  depth,  (c)  Vertical  and  temporal  distribution  for  a  pulse  500 
m  outbound  of  PS  (the  VRA  is  at  the  same  location  but  PS  is  500  m  closer  to  the  SRA). 

representation  of  the  field  produced  by  the  TRM.  Using  Eq.  33, 

(38) 

This  expression  can  be  used  to  show  that  the  TRM  produces  focusing  in  time  as  well  as  in  space. 
Focusing  in  time  occurs  because  a  form  of  matched  filtering  occurs.  To  understand  this,  examine  the 
TRM  field  at  the  focus  point  (that  is,  take  r  —  R,  z  =  Zpg  in  Eq.  38).  Neglecting  density  gradients, 
reciprocity  allows  the  interchange  Ga}{R)  Zps i  ^j)  =  Gu{RiZjyZps).  Then  the  time-domain  equivalent  of 
Eq.  38  is 

Ppcir,  t)  =  /  ^  [/  ^ 

where  the  time-domain  representations  of  the  Green’s  function  and  probe  pulse  are  used.  Note  that  the 
Green’s  function  is  correlated  with  itself.  This  operation  is  matched  filtering,  with  the  filter  matched 
to  the  impulse  response  for  propagation  from  the  probe  source  to  the  jih  SRA  element.  This  operation 
gives  focusing  in  the  time-domain,  that  is,  it  reduces  the  time  elongation  due  to  multipath  propagation 
[17],  The  sum  over  array  elements  is  a  form  of  spatial  matched  filtering,  analogous  to  that  employed  in 
the  matched  field  processor.  In  addition,  this  sum  further  improves  temporal  focusing  as  the  temporal 
sidelobes  of  the  matched  filters  for  each  channel  tend  to  average  to  zero  which  also  is  analogous  to 
broadband  matched-field  processing  results  [23],  Finally,  note  that  the  integral  over  in  Eq.  39  is 
a  convolution  of  each  matched-filtered  channel  impulse  response  with  the  time  reversed  and  delayed 
probe  pulse.  As  a  consequence,  this  pulse  is  not  matched  filtered,  for  example,  a  linear  FM  up-sweep 
will  appear  as  a  down-sweep  at  the  focus  and  will  not  be  compressed. 

Figure  18(a)  shows  a  simulation  for  a  50  ms  rectangular  pulse  with  center  frequency  445  Hz  for  the 
same  geometry  used  in  Fig.  17(a)  as  received  at  the  SRA  and  Fig.  18(b)  shows  the  pulse  as  transmitted 
to  a  plane  at  a  range  of  6.3  km,  the  range  of  PS.  Four  sources  were  excluded  from  the  simulation  because 
these  phones  were  not  used  in  the  experiment.  Note  the  temporal  focusing;  that  is,  the  50  ms  pulse 
disperses  to  about  75  ms  at  the  SRA  but  the  time  reversed  pulse  received  at  the  VRA  is  compressed 


^ pc(t' ,  i) 


i=i'' 
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Figure  19:  Experimental  results  for  probe  source  PS  and  VRA  at  same  range,  (a)  The  pulse  data 
received  on  the  SRA  for  PS  at  depth  of  40  m.  (b)  The  pulse  data  received  on  the  SRA  for  PS  at  depth 
of  75  m.  (c)  The  data  received  on  the  VRA  from  the  time  reversed  transmission  of  pulses  shown  in  (a). 
The  VRA  is  40  m  inbound  from  the  focus  as  determined  by  DGPS.  (d)  The  data  received  on  the  VRA 
from  the  time  reversed  transmission  of  pulse  shown  in  (b) .  The  VRA  is  40  m  outbound  from  the  focus 
as  determined  by  DGPS. 

(focused)  to  50  ms  as  opposed  to  exhibiting  even  further  time  dispersion.  On  the  other  hand  Fig.  18(c) 
shows  a  pulse  500  m  outbound  of  PS  (i.e.,  the  VRA  is  at  the  same  location  but  PS  is  500  m  closer  to 
the  SRA).  The  pulse  is  not  spatially  focused  and  it  is  temporally  more  diffuse  than  the  result  for  the 
focal  spot. 

3.2.3  Properties  of  the  Focal  Region 

The  TRM  focus  is  robust,  provided  the  SRA  adequately  samples  the  field  in  the  water  column.  First, 
the  focus  tends  to  depend  primarily  on  the  properties  of  the  ocean  near  the  focus  and  tends  to  be 
independent  of  (the  possibly  range- dependent)  properties  of  the  medium  between  the  SRA  and  the 
focus.  Temporal  changes  in  the  medium  due  to,  for  example,  surface  waves  and  internal  waves  degrade 
the  focus,  but  this  degradation  will  be  tolerable  if  the  average  (or  coherent)  Greens  function  is  not 
severely  reduced  by  these  time  variations.  Generally,  the  shape  of  the  focus  is  approximated  by  the 
field  that  a  point  source  placed  at  the  focus  generates  after  non-propagating  modes  are  subtracted. 
Thus,  if  absorption  or  scattering  tends  to  eliminate  high-order  modes,  the  focus  will  be  comprised  of 
the  remaining  lower  order  modes  and  will  be  relatively  broader.  Very  roughly,  the  vertical  width  of  the 
focus  will  be  equal  to  the  water  depth  (or  depth  of  the  duct)  divided  by  the  number  of  contributing 
modes  if  the  sound  speed  (in  the  duct)  is  not  strongly  dependent  upon  depth. 

The  TRM  focus  is  also  robust  with  respect  to  array  shape  [13]  provided  the  shape  does  not  change 
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Figure  20:  Experimental  results  for  50-ms,  445  Hz  center  frequency  probe  source  (PS)  at  81  m  depth 
and  various  ranges,  R,  between  the  PS  and  the  source-receive  array  (SRA).  (a)  i?=4.5  km.  (b)  R=7.7 
km.  (c)  il=15  km.  (d)  i?=20  km.  (e)  R=SQ  km.  (f)  R=15  km.  The  vertical  receive  array  (VRA)  is 
within  60  m  of  the  PS.  Both  the  VRA  and  PS  were  suspended  from  the  ALLIANCE  except  (f)  which 
was  from  a  rf-telemetered  VRA  at  15  km  range  in  other  run. 
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Figure  21:  The  energy  over  a  0.3-s  time  window  of  the  pulse  received  on  VRA  as  a  function  of  depth 
for  various  ranges.  The  depth  of  the  probe  source  was  81  m.  (a)  Experimental  results,  (b)  Simulation 
results.  Note  the  sidelobe  at  around  40  m  depth  and  15  km  range  in  both  cases. 

between  the  probe  reception  and  time-reversed  transmission.  This  property  makes  it  unnecessary  to 
know  the  exact  shape  of  the  TRM  array  and  offers  a  considerable  advantage  over  conventional  beam- 
forming.  The  focal  properties  are  discussed  in  detail  in  Sect  4. 

3.2.4  Variable  Range  Focusing 

The  time-reversal  mirror  technique  has  been  extended  to  refocus  at  ranges  other  than  that  of  the  probe 
source  [20].  The  technique  involves  retransmitting  the  data  at  a  shifted  frequency  according  to  the 
desired  change  in  focal  range,  such  that 

(Aw/w)  =  ^(AR/R),  (40) 

where  the  invariant  j3  is  determined  by  the  properties  of  the  medium.  Typically,  0  is  approximately 
equal  to  1  in  a  large  class  of  acoustic  waveguides,  which  means  that  a  5%  increase  in  frequency  results  in 
a  5%  increase  in  focal  range.  On  the  other  hand,  j3  fv -3  for  a  surface  sound  channel  or  in  a  waveguide 
with  a  deeply  submerged  channel  axis.  In  this  case,  retransmitting  frequency  should  be  decreased  to 
increase  the  focal  range.  The  frequency  shift  can  be  implemented  easily  in  near  real-time  by  a  FFT 
bin  shift  prior  to  retransmission.  Recently,  variable  focusing  in  range  has  been  applied  to  medical 
ultrasonics  in  free  space  [22]. 

3.3  Implementation  of  a  Time  Reversal  Mirror  in  the  Ocean 

Two  experiments  [19,  21]  were  conducted  by  the  SACLANT  Undersea  Research  Centre  (Dr.  Tuncay 
Akal  was  the  overall  Chief  Scientist  of  the  experiments)  and  the  Marine  Physical  Laboratory  of  Scripps. 
An  assortment  of  runs  was  made  to  examine  the  structure  of  the  focal  point  region  and  the  temporal 
stability  of  the  process.  The  TRM  process  was  successfully  demonstrated  out  to  a  range  of  thirty  km 
in  these  experiments.  Here  we  will  be  reporting  on  some  of  the  results,  (note  that  range  refers  to  the 
distance  from  the  SRA): 


Kuperman  and  Jackson:  Ocean  Acoustics 


27 


3,3.1  Demonstration  of  TRM  in  the  Ocean 

In  the  first  experiment  he  vertical  receive  array  VRA  was  deployed  at  a  range,  determined  by  DGPS, 
of  6.24  km  from  the  SRA  and  the  probe  source  PS  was  deployed  at  two  different  depths,  40  m  and  75 
m.  Figure  19  shows  the  pulse  as  received  on  the  SRA  and  VRA  for  both  source  depths.  The  data  at 
the  SRA  is  a  combination  of  signal  and  noise.  A  233  ms  window  was  digitized  and  time  reversed  for 
transmission  to  the  VRA.  When  the  VRA  and  PS  have  the  same  range  (experimentally  within  40  m  by 
a  DGPS  measurement)  to  the  SRA,  we  see  the  focusing  as  predicted  in  Section  3.1  for  a  probe  source 
at  40  m  depth  and  similar  results  for  a  probe  source  at  75  m  depth.  Clearly,  we  have  implemented  a 
time  reversal  mirror  focusing  at  the  range  and  depth  of  the  probe  source. 

In  the  second  experiment,  the  probe  source  was  deployed  out  to  a  range  of  30  km.  Figure  20  shows 
the  results  for  PS  at  an  81-m  depth  for  five  different  ranges  from  the  SRA:  4.5,  7.7,  15,  20,  and  30 
km.  As  expected,  the  temporal  focus  remain  compact  while  the  spatial  focus  broadens  with  range  due 
to  mode  stripping.  The  latter  can  be  seen  more  quantitatively  in  Fig.  21(a).  Simulations,  shown  in 
Fig.  21(b),  using  the  measured  environment  confirm  that  the  focal  structure  is  consistent  with  the  above 

theory. 

To  summarize,  the  TRM  is  still  quite  effective  at  30  km,  particularly  for  a  deep  probe  source  in 
the  stable  part  of  the  water  column.  The  focal  region  has  a  vertical  extent  of  about  25  m  (3  dB  down 
from  the  peak)  at  a  range  of  30-km,  which  corresponds  to  20%  of  the  waveguide  thickness  at  a  range  of 
about  250  waveguide  depths.  Although  we  did  not  map  out  the  radial  extent  of  the  focus,  a  simulation 
using  the  measured  environment  which  has  successfully  described  this  process  indicates  that  the  radial 
size  of  the  focus  was  on  the  order  of  800  m  at  a  30  km  range  and  300  m  at  a  5  km  range. 

3.3.2  Stability  Measurements 

Here  we  just  show  examples  of  stability  measurements  from  the  first  experiment.  Two  stability  data 
collection  periods  for  the  probe  source  depths  of  40  m  and  75  m  were  made  for  one  hour  and  two  hours, 
respectively  (the  lengths  of  the  runs  were  dictated  by  experimental  circumstance).  Figure  22  shows  the 
results  of  these  runs.  These  plots  indicate  that  the  focus  was  considerably  more  stable  for  the  deep 
probe  source  vs  the  shallower  probe  source  and  that  the  focus  is  broader  for  the  shallower  probe  source. 
This  is  consistent  with  the  measured  sound  speed  structure.  A  discussion  on  the  stability  of  the  TRM 
will  be  given  in  Sec.  5. 

3.3.3  TRM  with  variable  range  focusing 

Here  we  present  results  which  experimentally  confirm  a  technique  to  change  the  range  focus  of  a  TRM 
as  described  in  Sec.  3.3.3  [20].  During  the  experiment,  the  frequency  shift  was  implemented  in  near  real¬ 
time  by  a  simple  FFT  bin  shift  of  the  probe  source  data  received  by  the  SRA  prior  to  retransmission. 
Fig.  23(a)  shows  the  considerable  defocusing  for  the  PS  at  a  depth  of  68  m  when  the  VRA  was  600-m 
inbound  of  the  probe  source.  Fig.  23(b)-(d)  show  frequency  shifts  of  20  Hz,  25  Hz  and  32  Hz  with  the 
best  focus  resulting  from  the  32  Hz  shift.  It  is  interesting  to  note  that  this  latter  result  corresponds  to 
/9  =  1.4,  and  in  essence,  this  procedure  is  a  way  of  determining  /9. 

To  summarize,  it  is  possible  to  shift  the  focal  range  on  the  order  of  10%  of  the  nominal  range  of 
the  probe  source.  The  theory  on  which  this  shift  is  dependent  is  valid  only  over  a  frequency  range 
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Figure  22:  Results  on  stability  of  the  focal  region,  (a)  Pulse  arrival  structure  at  VRA  for  probe  source 
at  40-m  depth  averaged  over  1  h.  (b)  Pulse  arrival  structure  at  VRA  for  probe  source  at  75-m  depth 
averaged  over  2  h,  (c)  Mean  and  standard  deviation  of  energy  in  a  0.3-s  window  for  40-m  probe  source, 
(d)  Mean  and  standard  deviation  of  energy  in  a  0.3-s  window  for  75-m  probe  source. 


in  which  the  mode  shapes  do  not  change  significantly.  Frequency  shifts  of  greater  than  about  10% 
violate  this  condition.  A  practical  limitation  also  comes  from  the  transducer  characteristics  of  the  SRA 
whose  resonance  around  445  Hz  with  a  3  dB  bandwidth  of  approximately  35  Hz  as  shown  in  Fig.  B1  in 
Ref. [19].  Therefore,  it  is  difficult  to  excite  the  pulse  at  a  carrier  frequency  more  than  10%  offset  from 
the  original  resonance  frequency. 

3.4  Summary  of  TRM  experiments 

The  description  above  has  been  limited  mostly  to  the  implementation  of  a  time  reversal  mirror  in  the 
ocean.  The  experiments  also  studied  other  aspects  of  the  TRM  process  such  as  stability,  fluctuations, 
range  dependence,  out  of  plane  defocusing.  These  topics  will  be  discussed  elsewhere  in  this  workshop. 


4  The  range-dependent  ocean  waveguide 

The  factors  that  control  phase- conjugate  focusing  in  severely  range-dependent  environments  will  be 
examined  by  considering  a  general  non-uniform,  non-adiabatic  waveguide.  The  conditions  for  “ideal” 
phase-conjugate  focusing  in  such  a  waveguide  will  be  derived  and  this  will  implicitly  identify  the  factors 
that  degrade  focusing.  To  simplify  the  discussion,  only  vertical  source-receiver  arrays  will  be  considered. 
The  main  objective  is  to  generalize  Eq.  35  to  the  range-dependent  case,  using  the  approach  given  by 
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Figure  23:  Experimental  results  for  the  PS  at  a  depth  of  68  m.  (a)  Out-of-focus  results  on  the  VRA 
when  the  VRA  is  600  m  inbound  of  the  PS.  (b)  -20  Hz  frequency  shift,  (c)  -25  Hz  frequency  shift,  (d) 
-32  Hz  frequency  shift.  The  -32  Hz  shift  in  (d)  shows  the  best  focus,  which  corresponds  to  /S  =  1.4. 


Siderius,  et  al.  [32]  in  connection  with  the  “guide  source”  concept.  In  this  approach,  small  regions 
near  the  probe  source  and  SRA  are  assumed  to  be  range-independent,  but  the  larger  region  between  is 
allowed  to  have  arbitrary  range  dependence  in  bathymetry  and  sound  speed.  Losses  are  neglected  and 
will  be  discussed  later  in  qualitative  terms. 

The  Green’s  function  for  the  probe  field  near  the  probe  source  is  approximated  using  range-independent 
normal  modes.  /  \  /  d  \ 


^  ,/k„{R)\r-R\ 

Similarly,  the  Green’s  function  for  the  probe  field  at  the  SRA  is  written  in  the  form 

G.(0,  ;  fl, »..)  =  E 

„  y/Kn{^)K 

The  modal  eigenfunctions  in  the  vicinity  of  the  probe  source  and  SRA  are  denoted  Un{R.z)  and  tin(0,  r), 
respectively.  The  corresponding  eigenvalues  are  kn{R)  and  kn{0).  These  Greens  functions  do  not  bear 
the  superscripts  1  and  2  introduced  earlier  because  a  time-invariant  environment  is  under  consideration. 
The  subscript  u  is  used  here  in  the  same  sense  as  in  the  main  text.  The  mode  amplitudes  for  the  near- 
source  Green’s  function  are 

^g--iir/4 

071(2:^5)—  t  7^n(.R)  ^ps)  } 

y/S7rp{Zp,) 
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Figu're  24:  Coordinate  definitions  for  source-receiver  array  (SRA),  probe  source,  and  field  point  in 
range-dependent  ocean. 


and  the  mode  amplitudes  for  the  Green’s  function  near  the  SRA  are  given  by  the  linear  transformation 

^  .  (44) 

n 

For  convenience,  it  is  assumed  that  there  are  the  same  number  of  modes  near  the  source  and  near  the 
SRA,  so  that  Umn  is  a  square  matrix.  Cases  for  which  these  numbers  are  similar  but  not  equal  can  be 
treated  by  discarding  high-order  modes.  The  matrix  Umn  includes  any  mode  coupling  that  is  due  to 
the  range  dependence  of  the  ocean  and  is  defined  in  such  a  way  as  to  be  independent  of  source  depth. 
Furthermore,  to  the  extent  that  absorption  loss  in  the  water  column  and  seafloor  can  be  neglected,  Umn 
is  unitary: 

The  field  produced  by  the  SRA  is 

J 

<2:,  w)  =  0,  5  i?,  Zpj)  .  (do) 

i=i 


The  Green’s  function  for  propagation  from  the  jth  array  element  to  the  field  point  (r,  z)  can  be 
expressed  in  terms  of  the  Green’s  function  for  propagation  in  the  opposite  direction  by  using  reciprocity: 


P\^j  ) 


In  terms  of  mode  amplitudes, 


GUr,z-,0,  Zj)  =  4^:  y 
'  P(^i)  V 


(46) 


(47) 


where  the  mode  amplitudes,  c„(z),  are 


n» 


(48) 
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the  mode  amplitudes,  Cn(z),  are  essentially  the  same  as  the  bn{z)f  but  with  the  source  range 
i  shifted  by  r  —  i2. 

iOns  42  and  47  can  be  inserted  in  Eq.  45  to  obtain  an  expression  for  the  phase-conjugate  field 
-dependent  waveguide. 


D  / _ 


_  Um(0,2j)«n(0.^i)  (rr,\ 

-  '  (50) 

al  case,  the  SRA  spans  the  entire  water  column  with  elements  having  uniform  spacing,  c?a, 
lodal  eigenfunctions  have  negligible  amplitude  in  the  bottom.  In  this  case,  the  sum  over  array 
in  Eq.  50  approximates  the  orthonormality  integral  for  modal  eigenfunctions  (Eq.  17),  and 
in  be  taken  equal  to  S^n  •  This  ideal  can  be  approached  quite  closely  in  the  environment  of  the 
iiean  experiments.  Using  the  environmental  parameters  defined  in  Fig.  17,  and  considering 
irst  12  modes,  an  array  with  36  elements  with  spacing  da  =  3.33  m  and  with  the  shallowest 
44  m  below  the  surface  gives  diagonal  elements  in  Amn  da  that  are  within  3  percent  of  unity 
agonal  elements  that  are  of  order  0.03  or  less. 

ling  to  the  derivation  of  the  conditions  for  ideal  phase-conjugate  focusing,  take  =  Smnida 
to  obtain 


PpciT'i  Z\ui)  ^  —  7=  ^  ^  QmnQm(-2^)Qn(^pg)^ 


ik^(R){T-R) 


due  to  absorption  and  scattering  are  detrimental  to  phase-conjugate  focusing,  as  they  cause 
>n  of  higher-order  modes,  yielding  a  blurrier  focus  than  would  be  possible  with  lower  loss. 
)re,  this  blurring  will  increase  as  the  range  between  the  source  and  the  array  increases  owing 
ong  range-  and  mode-number  dependence  of  attenuation.  Figure  25  shows  the  loss  of  focus 
in  a  high-loss  silty  environment  as  compared  to  a  lower-loss  sandy  environment.  As  explained 
1.2,  losses  in  the  ocean  waveguide  are  primarily  due  to  refraction  of  energy  into  the  seafloor, 
ifloors  with  relatively  high  sound  speed  (sand)  are  less  lossy  than  seafloors  with  relatively  low 
ed  (silt).  As  a  result,  the  sandy  environment  supports  a  greater  number  of  acoustic  modes, 
a  sharper  focus. 

'loping  the  ideal  phase  conjugation  case,  losses  are  set  to  zero  and  the  mode  coupling  matrix, 
ken  to  be  unitary.  If  the  mode  dependence  of  ki(Q)  in  Eq.  52  is  neglected, 


hase-conjugate  field  for  an  ideal  array  in  a  lossless  environment  can  be  approximated  as 


Ppe(r,2;w)  =  ^ 


Sv  p{^Zpf)km(0^day/Rr 


Kuperman  and  Jackson:  Ocean  Acoustics 


32 


Focusing  Wth  Sand  Saafloor,  Run  rll 


Figure  25:  Comparison  of  focusing  for  sand  and  silt  seafloors.  The  sand  seafloor  has  sediment /water 
density  ratio  2.0  and  sound  speed  ratio  1.13,  with  absorption  of  0.135  dB  /  wavelength.  The  corre¬ 
sponding  numbers  for  the  silt  seafloor  are  1.5,  1.027  and  0.1.  The  gray  scale  is  linear  in  the  relative 
field  magnitude. 


‘Apart  from  inessential  factors,  this  expression  is  the  same  as  Eq.  35  which  was  derived  for  the  range- 
independent  case.  Even  though  Eq.  54  represents  the  ideal  case,  it  illustrates  properties  that  actual 
.phase-conjugate  arrays  may  possess,  provided  they  are  not  too  far  from  ideal.  One  such  property  is 
independence  of  the  focus  pattern  upon  the  distance  between  the  probe  source  and  the  array  (when 
absorption  can  be  neglected  and  apart  from  the  cylindrical  spreading  factor  l/\/^).  Even  more 
strikingly,  the  focus  field  is  independent  of  the  (possibly  range-dependent)  environment  between  the 
focus  and  the  array  (see  examples  presented  by  Siderius  et  al.  [32]  ).  That  is,  the  focus  depends 
only  on  the  local  properties  of  the  water  column  and  sea  floor  and  is  not  affected  by  bathymetry  or 
range-dependent  water-column  properties  in  the  region  between  the  array  and  the  focus,  provided  the 
latter  do  not  change  appreciably  during  the  two  propagation  cycles.  This  means  that,  in  the  ideal  case, 
phase  conjugation  is  not  affected  by  time-invariant  forward  scattering  due  to  bathymetry,  fronts,  etc. 
It  also  implies  that,  in  simulations  of  phase-conjugate  focusing,  it  is  important  to  accurately  model 
the  ocean  in  the  vicinity  of  the  focus,  but  less  accuracy  is  required  for  the  more  distant  parts  of  the 
propagation  path.  One  important  reservation  must  be  added  at  this  point.  The  derivation  above  is 
essentially  two-dimensional  in  that  cross-range  spatial  variation  of  the  ocean  is  neglected. 

Figure  26  shows  the  results  of  a  simulation  which  illustrates  the  invariance  of  the  focus.  In  this  figure, 
range-dependent  and  range-independent  environments  having  the  same  seafloor  type  (sand)  yield  very 
similar  focal  fields,  even  though  the  bathymetry  between  the  probe  source  and  SRA  is  greatly  different 
in  the  two  cases.  This  is  a  non-ideal  case  with  realistic  parameters,  including  loss.  Figure  27  shows  a 
case  in  which  losses  destroy  focus  invariance.  In  this  case,  the  probe  source  field  propagates  upslope 
to  reach  the  SRA.  Modes  are  lost  during  upslope  propagation,  leading  to  a  broader  focus  than  in  the 
range-independent  situation.  In  contrast,  modes  are  not  lost  in  the  downslope  propagation  of  Fig.  26. 

The  invariance  seen  in  the  ideal  ocean  waveguide  case  is  similar  to  that  predicted  for  the  ideal,  closed 
array,  except  that  the  time- reversed  field  is  a  standing  wave  in  the  latter  case.  In  the  present  case,  the 
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Figure  26:  Comparison  of  focusing  for  range-independent  and  range-dependent  (downslope  propagation 
of  probe  source)  environments.  The  sediment  acoustic  properties  are  the  same  as  those  of  the  sand 
seafloor  of  Fig.  25.  The  gray  scale  is  linear  in  the  relative  field  magnitude. 


probe  source  field  (including  only  propagating  modes)  is  given  by  Eq.  41  which  can  be  put  in  the  form 


G(r ,z\R,  ^pj)  — 


le 


-iirfA 


=E 


Un{R,  z)Un{R,  2:p,)e**’‘(^^l’'~^l 


(55) 


p{zp>)\/^A'r  -  R\  „  y/kniR) 

Apart  from  a  difference  in  spreading  loss  and  an  overall  phase  difference,  Eqs.  54  and  55  are  quite 
similar.  There  is  a  slight  term-by-term  difference  owing  to  differing  factors  involving  modal  eigenvalues, 
but  the  primary  difference  is  in  the  propagation  phase  factor.  The  source  field  propagates  ttway  from 
the  source  location  while  the  phase  conjugate  field  propagates  past  the  source  location  in  the  direction 
away  from  the  SRA. 


5  The  Effect  of  Ocean  Fluctuations  on  Phase  Conjugation 

The  ocean  is  random  in  several  respects.  The  surface  and  bottom  boundaries  are  randomly  rough, 
and  the  seawater  is  inhomogeneous  owing  to  stratification,  turbulence,  and  internal  waves.  In  addition, 
wave-generated  bubbles  introduce  inhomogeneity  in  a  layer  near  the  surface,  and  biological,  hydrody¬ 
namic,  and  geophysical  proces.ses  make  the  ocean  sediment  inhomogeneous.  Sound  scattering  due  to 
randomness  of  the  seafloor  is  time  independent,  and  is  not  expected  to  have  deleterious  effects  on  phase- 
conjugate  focusing.  In  fact,  the  range-independent  environment  discussed  earlier  can  be  considered  a 
kind  of  scattering  example,  and  it  was  seen  that  phase  conjugation  adapts  nicely  to  the  complexities 
of  this  case.  Generally,  phase  conjugation  can  be  expected  to  operate  robustly  in  the  presence  of 
time-independent  scattering.  Problems  arise  when  the  scattering  is  time-dependent,  since  the  medium 
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Figure  27;  Comparison  of  focusing  for  range-independent  and  range-dependent  (upslope  propagation 
of  probe  source)  environments.  The  sediment  acoustic  properties  and  gray  scale  are  the  same  as  those 

of  Fig.  25. 


changes  between  the  probe  and  retransmission  cycles.  This  is  the  case  with  scattering  by  surface  waves 
and  internal  waves,  for  which  one  would  like  to  know  the  characteristic  stability  time  over  which  a  gi\en 
probe  will  yield  satisfactory  focusing. 

5.1  Time-independent  volume  scattering 

As  noted  above,  volume  inhomogeneities  in  the  ocean  are  inherently  time  dependent,  owing  to  the 
presence  of  internal  waves.  These  waves  have  rather  long  periods  (minutes  to  hours),  however,  so 
that  volume  scattering  can  be  approximated  as  time  independent  if  the  time  between  probing  and 
retransmission  is  sufficiently  short.  In  this  case  one  does  not  expect  scattering  to  degrade  the  focus, 
in  fact,  it  will  be  seen  that  the  focus  is  improved  by  scattering.  In  static  environments,  the  subscripts 
1  and  2  can  be  dropped  in  Eq.  12  because  there  is  no  change  in  the  environment  between  the  probe 
transmission  and  the  time-reversed  retransmission.  It  is  useful  to  consider  the  formal  average  of  the 
focus  field,  imagining  an  ensemble  of  time-reversal  measurements  made  using  a  hypothetical  ensemble 
of  random  environments.  Often  in  underwater  acoustics,  the  average  of  a  field  over  an  ensemble  of 
realizations  is  not  comparable  in  magnitude  to  the  field  seen  in  single  realizations,  because  scattering 
often  causes  fluctuations  that  are  comparable  to  or  greater  than  the  mean.  This  is  not  the  situation  with 
phase  conjugation  if  the  SRA  spans  an  appreciable  portion  of  the  water  column,  as  the  arguments  made 
in  connection  with  Eq.  55  show  that  the  focus  field  is  not  very  strongly  dependent  on  perturbations  of 
the  environment.  Thus,  it  is  useful  to  examine  the  average  focus  field,  taking  it  as  an  approximation  to 
the  field  seen  in  any  given  realization. 


Figure  28:  Coordinate  definitions  for  static  volume  scattering  environment. 

J 

<  Ppeir; u)  >=  ^  K(v,  Vp, , r^)  ,  (56) 

i=i 

where 

K{r,Tps,rj)  =<  Ga;(r,rj)G*(rp3,rj-)  >  .  (57) 

The  first  moment  of  the  focus  field  depends  on  the  second  moment  of  the  Greenes  function.  In  some 
simple  cases,  the  second  moment  of  the  Green’s  function  is  of  the  form  (see  [33]  for  an  example  and 
references,*  and  [34]  for  dynamic  case) 

j  0, Gojw(j’, -2:,  0,  2:j)(!j0  0. -Sj )e  ^  (^^) 

where 

H  =  ^J^''‘{l-p[AziIt!)]}dR'  ,  (59) 

and  where  Go,w(^,  0,  zj)  is  the  Green’s  function  for  the  average  medium,  />( 2:)  is  the  vertical  correlation 

function  for  refractive  index  (p(0)  =  1),  is  the  phase  variance  for  a  path  of  length  R,  and  Az{R')  is 
vertical  separation  of  ray  paths  (Fig.  5.1). 

Using  this  result,  the  mean  phase-conjugate  field  is 

<  Ppc{R,  -2:;a;)  >=  Ppco{R,  z]uj)e~‘^  (60) 

where 

j 

Ppcoi^Ri  ^ }  ^)  ^  ^  G^d,a'(^}  •2^?  0,  “2^^  )^0,a?(-^)  >  0,  ) 

i=i 


(61) 


Kuperman  and  Jackson:  Ocean  Acoustics 


36 


is  the  phase-conjugate  field  for  the  average  medium.  Following  Uscinski  and  Reeve  [35],  expand  the 
correlation  function  for  refractive  index  fluctuations  in  a  power  series: 

piz)  =  l-^+...  (62) 

*V 

Assuming  that  the  average  medium  is  infinite  and  non-refr active,  the  ray  paths  are  straight  lines  and 
the  integral  for  H  can  be  performed. 

«  =  ,e3) 

The  factor  e~'^  has  a  Gaussian  shape  in  z  centered  on  the  probe  source  with  1/e  full  width 


In  some  cases,  this  width  can  be  less  than  the  focal  spot  width  for  the  average  medium  [14],  [36].  This 
sharpening  of  the  focus  has  been  observed  experimentally  but  in  a  stronger  scattering  regime  than  the 
present  formalism  can  describe  [37]. 

5.2  Time-dependent  scattering  by  surface  waves 

Time-dependent  forward  scattering  due  to  surface  waves  causes  degradation  of  phase-conjugate  focusing. 
Following  Ref.  [19],  the  Green’s  function  will  be  decomposed  into  coherent  (unscattered)  and  incoherent 
(scattered)  parts: 

Ga(r,  r')  =  G(r,  r')  -f  ^Ga(r ,  r')  .  (65) 

The  subscript  a  takes  on  the  values  1  and  2  for  the  probe  and  conjugate  transmission  cycles,  respectively. 
The  coherent,  or  mean,  Green’s  function,  G(r,r')  is  not  assigned  a  subscript  because  the  random  time 
variations  axe  assumed  to  be  stationary  in  the  statistical  sense.  It  will  be  assumed  that  sufficient  time 
has  elapsed  between  the  probe  and  conjugate  transmission  cycles  that  variations  in  the  two  Green’s 
functions  are  uncorrelated. 

<  ^G2(rd,rc)(5Gt(r6,ra)  >  =  <  <5G2(rd,rc)6Gi(r6,ra)  >  =  0  .  (66) 

This  condition  was  very  likely  satisfied  in  the  Mediterranean  experiments  with  respect  to  scattering 
by  surface  waves,  which  have  correlation  time  scales  on  the  order  of  seconds,  while  the  time  between 
transmission  cycles  was  measured  in  minutes  and  hours. 

Combining  Eqs.  30,  65,  and  66,  the  mean  phase-conjugate  field  is 

j 

Ppc{r,  z-,u)=Y^  G(r,  Tj)G* {vj ,  rp,)  ,  (67) 

i=i 

and  the  variance  of  the  field  is 


lPpe(r,z;w)|2 _  |Pp^(r,z;w)p  =  ^  ^[G(r,rj')^ (r,rj/)A"jj-(rp,)+ 
G{Tps ,  rj)G  (xp, , Yji)Kjji{r)  +  » 


Kuperman  and  Jackson:  Ocean  Acoustics 


37 


Figure  29:  Vertical  slice  through  mean  focus  for  selected  RMS  surface  wave  heights. 


where 

J^jj‘(r)=<SGa(rj,r)SGa(rj>,r)>  .  (69) 

The  covariance,  Kjjj(r)y  is  proportional  to  the  correlation  between  the  incoherent  field  at  elements  j 
and  f  of  the  array,  with  a  unit  point  source  situated  at  r.  In  deriving  Eq.  68,  free  use  was  made  of 
reciprocity  (which  allows  interchange  of  the  two  arguments  of  the  Green’s  function)  and  stationarity 
(which  means  that  SGi  and  SG2  have  identical  statistics). 

Equations  67  and  68  are  general  and  include  three-dimensional  scattering  (i.e.,  in-plane  and  out- 
of-plane  scattering).  They  lead  to  two  general  conclusions  regarding  focusing  in  ocean  experiments 
for  those  cases  in  which  sufficient  time  elapsed  between  the  two  transmission  cycles.  First,  the  mean 
focus  field,  that  is,  the  focus  field  averaged  over  many  independent  probe-conjugate- transmission  cycles, 
is  obtained  by  using  the  coherent  Green’s  function  in  place  of  the  actual  (random)  Green’s  function. 
Second,  and  most  important,  the  field  near  the  focus  does  not  fluctuate  appreciably,  that  is,  it  is  well 
approximated  by  the  mean  focus  field.  As  shown  in  [19],  careful  inspection  of  Eq.  68  shows  that  the 
variance  of  the  phase-conjugate  field  is  not  localized  near  the  focus,  but  is  spread  diffusely  in  range 
and  depth.  Thus,  near  the  focus,  the  mean  field  dominates,  unless  scattering  is  strong  enough  to 
diminish  the  mean  Green’s  function  to  such  a  degree  that  focusing  is  essentially  destroyed.  It  is  rather 
straightforward  to  compute  the  mean  field  in  the  presence  of  surface-wave  scattering  [38].  Figure  29  is 
an  example  of  such  a  calculation  using  the  modal  propagation  code,  KRAKEN  [39]  with  parameters 
appropriate  to  the  1996  Mediterranean  experiment.  For  RMS  wave  heights  within  the  range  measured 
during  the  experiment,  the  degradation  of  the  mean  focus  is  predicted  to  be  slight.  This  agrees  with  the 
experimental  fact  that  the  measured  focus  was  stable  on  the  short  time  scales  characteristic  of  surface 
wave  motion. 

5.3  Time-dependent  scattering  by  internal  waves 

Internal  waves  are  ubiquitous  in  the  ocean  and  arise  because  of  density  stratification.  Like  surface  waves, 
internal  waves  involve  an  interplay  of  kinetic  energy  and  gravitational  potential  energy.  Internal  waves 
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cause  a  slow  orbital  motion  of  the  water  which  distorts  the  sound  speed  profile  in  both  the  vertical  and 
horizontal,  leading  to  random  inhomogeneity  in  the  acoustic  index  of  refraction.  This  inhomogeneity, 
in  turn,  causes  time-dependent  fluctuations  in  sound  propagation,  which  affect  the  stability  of  phase 
conjugation. 

In  this  section,  a  brief  summary  of  relevant  internal  wave  properties  is  given,  concentrating  on  the 
statistical  description  needed  for  modeling  of  propagation  fluctuations.  For  more  complete  accounts 
of  internal  wave  dynamics  and  statistics,  see  [40]  and  [41].  For  treatments  of  the  problem  of  acoustic 
scattering  by  internal  waves,  see  [42,  43,  44]. 

The  vertical  displacement  due  to  internal  waves,  C(r,  t),  can  be  written  as  a  sum  over  normal  modes 

C(r,0  =  (70) 

K  i=i 

where  r  =  (R, 2:),  R  =  {x,y),  K  =  {Kx,Ky).  The  frequencies,  fi(K),  and  eigenfunctions,  Wj{z),  are 
obtained  by  solving  the  eigenvalue  equation 

(//  ”  +  K^[N\z)  -  fJ]W^  =  0  (71) 

subject  to  the  conditions  that  the  displacement  field  vanishes  at  the  surface  {z  =  0)  and  bottom  (z  =  A). 
The  parameter  //  is  the  inertial  frequency: 

//  =  2sm(^x)/T  (72) 

Here,  is  the  latitude  and  T  is  the  rotation  period  of  the  earth. 

The  function  N{z)  is  the  depth-dependent  buoyancy  frequency,  determined  by  the  vertical  gradient 
of  density. 

N\z)  =  _  ^^(,)]  (73) 

where  Pa{z)  is  the  density  profile  in  an  adiabatically  mixed  ocean,  and  g  is  the  acceleration  of  gravity. 
The  eigenfunctions  obey  the  following  orthonormality  relation: 

t  [[iv'(2)  -  mWj {z)Wp{z)dz  =  (74) 

-^0 

Figure  30  shows  a  buoyancy  frequency  profile  obtained  from  oceanographic  data  in  the  1997  Mediter¬ 
ranean  experiment  [45].  Several  of  the  resulting  eigenfunctions,  computed  for  a  particular  wavenumber 
are  also  shown  in  the  figure. 

The  space-time  covariance  of  the  displacement  field  is  defined  as  the  following  average, 

<  C(Ro  -f-R, +  r)C(Ro,  Z2,t)  >=  K{R,zi,Z2,  r)  (75) 

where  it  is  assumed  that  the  internal  wave  field  is  statistically  isotropic  and  stationary  in  the  transverse 
directions  and  stationary  with  respect  to  time  so  that  the  covariance  depends  on  /^  =  y/x^  -f  and  r, 
where  (a;,  y)  is  the  horizontal  vector  separation  of  the  two  points  for  which  the  covariance  is  measured. 
The  displacement  field  is  not  stationary  in  the  vertical,  hence  the  covariance  depends  on  both  zy  and 
Z2-  A  useful  general  expression  for  the  covariance  can  be  obtained  by  squaring  and  averaging  Eq.  70 
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Figure  30:  Buoyancy  frequency  profile  and  internal  wave  eigenfunctions  for  the  1997  Mediterranean 
experiment.  The  eigenfunctions  are  computed  for  horizontal  wavenumber  K  =  0.001  m”^. 


rOO 

KiR,  zi,Z2,  r)  =  V  /  Fif,j)Wjiz,)Wj{z2)MKR)cos{2wfT)df  (76) 

The  eigenfunctions,  Wj{z)^  and  the  horizontal  wavenumber,  depend  upon  both  the  mode  number, 
j,  and  the  frequency,  /. 

Henyey  et  al.  [42]  provide  a  shallow- water  form  of  Garrett-Munk  internal- wave  statistics  as  follows: 


(77) 


The  two  parameters  b  and  Egm  appear  as  a  product,  and  this  product  (having  units  of  length)  is 
treated  as  a  single  parameter  determining  the  energy  of  the  internal  wave  field.  Typically,  bEoM  is  of 
order  1  m.  The  parameter  n,  is  chosen  so  that  the  sum  of  H{j)  over  all  modes  is  unity.  The  parameter 
determines  which  modes  are  most  important.  In  deep  water,  =  3  is  typical,  but  smaller  values  are 
expected  in  shallow  water  [42]. 

Simulated  internal  wave  displacement  fields  can  be  generated  by  assigning  random  values  to  the 
amplitudes  ay(K)  in  Eq.  70.  The  resulting  random  sound  speed  field  follows  from  knowledge  of  the 
sound  speed  profile  at  zero  displacement,  which  may  be  estimated  using  an  average  over  several  measured 
profiles.  Figure  32  shows  simulations  of  focusing  using  a  probe  sent  at  t  =  0.  These  simulations,  for 
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Wavenumber  (1/m) 

Figure  31:  Internal  wave  frequency  vs.  wavenumber  for  the  1997  Mediterranean  experiment. 

an  acoustic  frequency  of  400  Hz,  employed  internal  wave  parameters  =  1  and  bEoM  =  2.0  m. 
This  value  of  bEoMi  the  internal  wave  energy  parameter,  is  much  larger  than  the  value  suggested  by 
oceanographic  data  gathered  during  the  1997  Mediterranean  experiment,  bEoM  =  0.135  m  [45].  The 
smaller,  more  realistic  value  produced  negligible  internal  wave  effects,  consistent  with  the  observed 
stability  of  the  focus.  With  the  stronger  internal  wave  field,  Fig.  32  shows  that  the  initial  sharp  focus 
degrades  substantially  on  time  scales  of  several  minutes. 

6  Conclusions 

The  ocean  is  an  acoustic  waveguide  whose  fundamental  properties  are  well  understood.  Prediction  of 
acoustic  phenomenon  requires  extensive  ocean  acoustic  environmental  knowledge  that  is  often  difficult  or 
impossible  to  obtain.  Matched  field  processing,  a  passive  localization  technique  is  dependent  on  accurate 
acoustic  information  to  construct  “replicas”  and  therefore  is  extremely  sensitive  to  mismatch.  Phase 
conjugation  and  time  reversal  mirrors  has  been  shown  to  be  the  active  implementation  of  a  matched 
field  processor  where  an  array  of  sources  is  used  to  propagate  the  measured  data  on  an  array  back  to  its 
origin.  In  essence,  the  ocean  is  used  to  construct  the  replica’s,  eliminating  the  possibility  of  mismatch  in 
a  reasonably  stable  ocean  environment.  Probably  the  most  amazing  aspect  of  implementing  a  TRM  in 
the  ocean  over  distances  as  great  as  30  km  is  that  the  performance  of  the  ocean  TRM  is  almost  exactly 
the  same  as  that  of  a  scaled  (MHz  regime)  ultrasonic  TRM  implemented  in  a  totally  stable  laboratory 
environment  out  to  distances  of  tens  of  centimeters  [24].  The  TRM  process  is  extremely  robust. 
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Appendix  A:  Parabolic  Equation  Model  (PE) 

The  PE  method  was  introduced  into  ocean  acoustics  and  made  viable  with  the  development  of  the  “split 
step  algorithm”  which  utilized  fast  Fourier  transforms  at  each  range  step  [25].  Subsequent  numerical 
developments  greatly  expanded  the  applicability  of  parabolic  equation. 


Standard  PE-Split  Step  Algorithm 

The  PE  method  is  presently  the  most  practical  and  encompassing  wave- theoretic  range-dependent 
propagation  model.  In  its  simplest  form,  it  is  a  far-held  narrow-angle  ±20°  with  respect  to  the 
horizontal — adequate  for  most  underwater  propagation  problems)  approximation  to  the  wave  equa¬ 
tion.  Assuming  azimuthal  symmetry  about  a  source,  we  express  the  solution  of  Eq.  10  in  cylindrical 
coordinates  in  a  source  free  region  in  the  form 

G(r,z)  =  -0(^2)  *  J(r),  (80) 

and  we  define  K^(r,  z)  =  KqTi^,  n  therefore  being  an  “index  of  refraction”  co/c,  where  cq  is  a  reference 
sound  speed.  Substituting  Eq.  80  into  Eq.  10  in  a  source  free  region  and  taking  Kq  as  the  separation 
constant  ,  J  and  0  satisfy  the  following  two  equations: 

d^J  IdJ  ^  .... 

+  =0;  (81) 
r  dr 

^ 

Equation  81  is  a  Bessel  equation  and  we  take  the  outgoing  solution,  a  Hankel  function,  H^iKor)^ 
in  its  asymptotic  form  and  substitute  it  into  Eq.  82,  together  with  the  “paraxial”  (narrow  angle) 
approximation, 

/OQ. 


to  obtain  the  parabolic  equation  (in  r), 


0  +  2i/<o|^  +  /’£:2(n"-l)i/-  =  O,  (84) 

where  we  note  that  n  is  a  function  of  range  and  depth.  We  use  a  marching  solution  to  solve  the 
parabolic  equation.  There  has  been  an  assortment  of  numerical  solutions  but  the  one  that  still  remains 
the  standard  is  the  so  called  “split  step  algorithm”  [25]. 

We  take  n  to  be  a  constant;  the  error  this  introduces  can  be  made  arbitrarily  small  by  the  appropriate 
numerical  gridding.  The  Fourier  transform  of  0  can  then  be  written  as 


1 

^  y  exp(-2S2r)  dz 


which  together  with  Eq.  84  gives 


-  s^x  ±  ±  -  1)X  =  0. 


The  solution  of  Eq.  86  is  simply 


x(r,s)  =  x(ro,5)  exp 


Kl{n^ 


-  1)  - 


(r-ro) 
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with  specified  initial  condition  at  tq.  The  inverse  transform  gives  the  field  as  a  function  of  depth, 

x(ro,s)exp  -^(n^  -  l)Ar  exp  exp(is.-)  ds,  (88) 

where  Ar  =  r—rQ.  Introducing  the  symbol  ^  for  the  Fourier  transform  operation  from  the  ^-domain  (as 
performed  in  Eq.  85)  and  as  the  inverse  transform,  Eq.  88  can  be  summarized  by  the  range-stepping 
algorithm, 

^(r  +  Ar,  z)  =  exp  -  l)Ar]  [(expC-^^^))  •  ^  z)^  ,  (89) 

which  is  often  referred  to  as  the  “split  step”  marching  solution  to  the  PE.  The  Fourier  transforms  are 
performed  using  FFT’s.  Equation  89  is  the  solution  for  n  constant,  but  the  error  introduced  when 
n  (profile  or  bathymetry)  varies  with  range  and  depth  can  be  made  arbitrarily  small  by  increasing 
the  transform  size  and  decreasing  the  range-step  size.  It  is  possible  to  modify  split  step  algorithm  to 
increase  its  accuracy  with  respect  to  higher  angle  propagation  [26]. 


Generalized  or  Higher  Order  PE  Methods 

Methods  of  solving  the  parabolic  equation,  including  extensions  to  higher  angle  propagation,  elastic 
media,  and  direct  time  domain  solutions  including  nonlinear  effects  have  recently  appeared  (see  Refs. 
[27,  28]  for  additional  references).  In  particular,  accurate  high  angle  solutions  are  important  when 
the  environment  supports  acoustic  paths  that  become  more  vertical  such  as  when  the  bottom  has  a 
very  high  speed  and  hence,  large  critical  angle  with  respect  to  the  horizontal.  In  addition,  for  elastic 
propagation,  the  compressional  and  shear  waves  span  a  wide  angle  interval.  Finally,  Fourier  synthesis 
for  pulse  modeling  requires  high  accurate  in  phase  and  the  high  angle  PE’s  are  more  accurate  in  phase, 
even  at  the  low  angles. 

Equation  82  with  the  second  order  range  derivative  which  was  neglected  because  of  Ineq.  83  can  be 
written  in  operator  notation  as 

[p2  -f  2iKoP  -b  K^iQ^  -  =  0,  (90) 


where 

Factoring  Eq.  90  assuming  weak  range  dependence  and  retaining  only  the  factor  associated  with  outgoing 
propagation  yields  a  one  way  equation 

Pxl>  =  iKo{Q  -  (92) 

which  is  a  generalization  of  the  parabolic  equation  beyond  the  narrow  angle  approximation  associated 
with  Ineq.  83.  If  we  define  Q  =  i/l  +  5  and  expand  Q  in  a  Taylor  series  as  a  function  of  q,  the  standard 
PE  method  is  recovered  by  Q  «  1  +  .69.  The  wide  angle  PE  to  arbitrary  accuracy  in  angle,  phase, 
etc.can  be  obtained  from  a  Pade  series  representation  of  the  Q  operator  [27] , 


Q = v^i + 9 = 1 +y^  Y 

i=i 


+  bj,nq 


+  0{q^ 


(93) 
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where  n  is  the  number  of  terms  in  the  Fade  expansion  and 


— 


2n  -|- 1 


JTT 


2?!  4“  1 


(94) 


The  solution  of  Eq.  92  using  Eq.  49  has  been  implemented  using  finite  difference  techniques  for  fluid 
and  elcLstic  media  [27].  A  “split  step”  Fade  algorithm  [29]  has  recently  been  developed 'which  greatly 
enhances  the  numerical  efficiency  of  this  method. 
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Appendix  B:  Units 

The  decibel  (dB)  is  the  dominant  unit  in  underwater  acoustics  and  denotes  a  ratio  of  intensities  (not 
pressures)  expressed  in  terms  of  a  logarithmic  (base  10)  scale.  Two  intensities,  /i  and  I2  have  a  ratio, 
I1/I2  in  decibels  of  lOlogi'i/72  dB.  Absolute  intensities  can  therefore  be  expressed  by  using  a  reference 
intensity.  The  presently  accepted  reference  intensity  is  based  on  a  reference  pressure  of  one  micropascal 
(pPa):  the  intensity  of  a  plane  wave  having  an  rms  pressure  equal  to  10“®  dynes  per  square  centimeter. 
Therefore,  taking  1  pPa  as  /2j  a  sound  wave  having  an  intensity,  of,  say,  one  million  times  that  of  a 
plane  wave  of  rms  pressure  1  jiPa  has  a  level  of  I01og(10®/1)  =  60  dB  re  1  pPa.  Pressure  (p)  ratios 
are  expressed  in  dB  re  1  pPa  by  taking  201ogpi/p2  where  it  is  understood  that  the  reference  originates 
from  the  intensity  of  a  plane  wave  of  pressure  equal  to  1  /iPa. 

The  average  intensity,  I,  of  a  plane  wave  with  rms  pressure  p  in  a  medium  of  density  p  and  sound 
speed  c  is  /  =  In  seawater,  pc  is  1.5  x  10®  so  that  a  plane  wave  of  rms  pressure  1 

dyne  j err?  has  an  intensity  of  0.67  x  10“^^PI^/cm^.  Substituting  the  value  of  a  micropascal  for  the  rms 
pressure  in  the  plane  wave  intensity  expression,  we  find  that  a  plane  wave  pressure  of  1  pPa  corresponds 
to  an  intensity  of  0.67  x  lO^^^kP/cm^  (i.e.,  0  dB  re  1  pPa). 
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